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Large Eddy Simulation and Direct 
Numerical Simulation of High Speed 
Turbulent Reacting Flows 


V. Adumitroaie* S.H. Frankelj C.K. Madnia and P. Givi 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, New York 14260-4400 


Abstract 

The objective of this research is to make use of Large Eddy Simulation (LES) and 
Direct Numerical Simulation (DNS) for the computational analyses of high speed re- 
acting flows. Our efforts in the first phase of this research conducted within the past 
three years have been directed in several issues pertaining to intricate physics of turbu- 
lent reacting flows. In our previous 5 semi-annual reports submitted to NASA LaRC, 
as well as several technical papers in archival journals, the results of our investigations 
have been fully described. In this progress report which is different in format as com- 
pared to our previous documents, we focus only on the issue of LES. The reason for 
doing so is that LES is the primary issue of interest to our Technical Monitor and that 
our other findings were needed to support the activities conducted under this prime 
issue. The outcomes of our related investigations, nevertheless, are included in the ap- 
pendices accompanying this report. The relevance of the materials in these appendices 
are, therefore, discussed only briefly within the body of the report. 

Here, results are presented of a priori and a posteriori analyses for validity assess- 
ments of assumed Probability Density Function (PDF) methods as potential subgrid 
scale (SGS) closures for LES of turbulent reacting flows. Simple non-premixed reacting 
systems involving an isothermal reaction of the type A + B Products under both 
chemical equilibrium and non-equilibrium conditions are considered. A priori analyses 
are conducted of a homogeneous box flow, and a spatially developing planar mixing 
layer to investigate the performance of the Pearson Family of PDF’s as SGS models. 


*One of the two Primary Graduate Research Assistants supported partially under this Grant. 

* Former Graduate RA supported partially by this Grant. Present Affiliation: Assistant Professor of 
Mechanical Engineering, Purdue University, West Lafayette, India 47109-1077. 
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A posteriori analyses are conducted of the mixing layer using a hybrid one-equa ion 
Smagorinsky /PDF SGS closure. The Smagorinsky closure augmented by the solution 
of the subjid turbulent kinetic energy (TKE) equation is employed to account for 
hydrodynamic fluctuations, and the PDF is employed for modeling ^eeffectsofscalar 
fluctuations. The implementation of the model requires the knowledge of the local 
values of the first two SGS moments. These are provided by additional modeled trans- 
"aits. In both « priori and u posteriori analyte., the predicted I 
appraised by comparison with subgrid averaged results generated by DNS. Baaed 
these results, the paths to be followed in future investigations are identified. 

TuTresearch is sponsored by the NASA Langley Resell. Center r unde. rGM 
NAG-1-1122. Dr. J. Philip Drummond, Theoretical Flow Physics Branch (TFP ), 
Mail Stop 156, Tel: 804-864-2298 is the Technical Monitor of this Grant 


1 Introduction 

In the last three decades computer simulation has gained a substantial impetus aa one of 
the most promising tools in scientific research, already having a major effect on the hu- 
man understanding and discovery of complex physical phenomena. Factually, in areas where 
experiments are very difficult, if not impossible, to be performed plasma physics, as- 

trophysics, cosmology, etc.), numerical simulation have become the only method of verifying 
theories. Turbulence, as one of the unsolved problems with significant range of applications, 
has been the subject of broad numerical investigations - the rate significantly increased 
with the advent of the supercomputer technology since early 80’s. In general, presently 
there are three methodologies by which turbulent flows are treated by computer simula- 
tions: Direct Numerical Simulations (DNS), Large Eddy Simulations (LES) and Reynolds 
Averaged Navier-Stokes simulations (RANS). Each of these portray its advantages and its 
(often severe) limitations. In a recent review,' DNS and LES have been grouped together 
as “model-free simulations", since both approaches are aimed to provide a detailed descrip- 
tion of the large structures of the flow evolution. Although DNS involves computation in 
the entire bandwidth associated with the turbulent field, it remains confined to relatively 
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low Reynolds numbers ■ to say the least - for the proximate future. Nevertheless, DNS has 

become indispensable in turbulence research in assessing and calibrating models and also in 
exploring new concepts . 2,3 

LES is considered somewhere between DNS and RANS computation.'-'* Over the past thirty 
years, since the early work of Smagorinsky ' 7 there has been relatively little effort, compared 
to that in RANS calculations, to make full use of LES for engineering applications. The 
most prominent model has been the Smagorinsky based eddy viscosity closure which relates 
the unknown subgrid scale (SGS) Reynolds stresses to the local large scale rate of strain. 
This viscosity is aimed to provide the role of mimicking the dissipative behavior of the 
unresolved small scales. With sound SGS models, LES is expected to be used extensively in 
the study of high Reynolds numbers turbulence and flow prediction in complex situations. 
The nonlinear trait of scales interaction makes the task of developing correct SGS models 
very challenging even in the case of non-reacting turbulent flows. The difficulty is augmented 
with the increase in the number of degrees of freedom variable density flows, and 
currents with finite electrical conductivity) and with coupling between different phenomena 
(such as fluid dynamics and chemistry in combusting flows).'* Dealing with this complexity 
necessitates serious efforts by the scientific community, consequently at this time the grand 
ensemble of turbulence theories tends to exhibit a somewhat familiar disorderly appearance. 

The extensions to one-equation models, typically based on the SGS turbulent kinetic energy, 18 ' *> 
have shown some improvements 2 ''* 1 over Smagorinsky based closures. This is particularly 
the case in simulating transitional flows where the assumption of a balance between the 
production and the dissipation of turbulent kinetic energy may not always be valid. Thus, 
the higher degree of freedom provided by one-equation closures allow more flexibility for the 
subgrid scale eddies to adjust to local flow conditions. 

A survey of combustion literature reveals relatively little work in LES of chemically reacting 
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turbulent flows . 1 - 23 It appears that Schumann 24 was one of the first to conduct LES of a 
reacting flow. However, the assumption made in this work to simply neglect the contribution 
of SGS scalar fluctuations to the filtered reaction rate is debatable. The importance of 
such fluctuations is well recognized in RANS of reacting flows in both combustion 25-27 and 
chemical engineering 28-31 problems. Therefore, it is natural to expect that these fluctuations 
will also have a significant influence in LES. 

Modeling of scalar fluctuations in RANS has been the subject of intense investigations since 
the pioneering work of Toor . 32 One approach which has proven particularly useful is based on 
the Probability Density Function (PDF) or joint PDF of scalar quantities . 33-38, 1,37,23 This 
approach offers the advantage that all the statistical information pertaining to the scalar 
field is embedded within the PDF. Therefore, once the PDF is known the effects of scalar 
fluctuations are easily determined. Because of their capabilities, PDF methods have been 
widely used in RANS for a variety of reacting systems (see Refs . 23 ’ 38 for recent reviews). A 
systematic approach for determining the PDF is by means of solving the transport equation 
governing its evolution . 39,40 In this equation the effects of chemical reactions appear in a 
closed form. However, modeling is needed to account for transport of the PDF in the do- 
main of the random variables. This transport describes the role of molecular action on the 
evolution of PDF. In addition, there is an extra dimensionality associated with the composi- 
tion domain which must be treated. These two problems have constituted a stumbling block 
in utilizing PDF methods in practical applications, and developments of turbulence closures 
and numerical schemes which can effectively deal with these predicaments have been the 
subject of numerous investigations within the past two decades. 

An alternative approach in PDF modeling is based on assumed methods. In these methods 
the PDF is not determined by solving a transport equation. Rather, its shape is assumed a 
priori usually in terms of the low order moments of the random variable(s). Obviously, this 
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method is ad hoc but it offers more flexibility them the first approach. Presently, the use 
of assumed methods in RANS is justified in cases where there is a strong evidence that the 
PDF assumes a particular distribution. 

Among these two approaches obviously the first one is preferable if an appropriate closure 
is available to account for the effects of molecular action. In its application in RANS, 
traditionally the family of models based on the coalescence/dispersion (C/D) closures, 41-43 
or least mean square methods 44 have been employed. These closures are plausible from 
a computational standpoint and can be effectively simulated via Monte Carlo numerical 
methods. 45 However, there are several drawbacks associated with these closures that restrict 
their use for accurate and reliable predictions. 46 " 48 Some of these drawbacks are overcome in 
the newly developed Amplitude Mapping Closure (AMC). 49,50,23 This has been established 
in a number of recent validation assessments of the AMC by means of comparison of its 
predicted results with those of DNS, 51-54 and experimental 55 data. 

Despite its demonstrated properties, there are some deficiencies associated with the AMC 
which require further investigations. These are discussed in detail in Ref.; 56 the most serious 
of these are: (1) the “single- point” nature of the closure, (2) the difficulties associated with 
its numerical implementation, and (3) its inability to account for the migration of scalar 
bounds as mixing proceeds. The first problem is shared with C/D models and indicates the 
deficiency of the approach in accounting for the variation of turbulent length (time) scales. 
The other problems are exclusive to AMC and can cause difficulties in its applications. 

Considering the current state of affairs in PDF modeling, it can be cautiously concluded 
that assumed PDF methods are somewhat more “feasible” than the transport equation ap- 
proach for simulating practiced problems. This is not to advocate the superiority of assumed 
methods. Rather, it is to encourage further research on the first approach before it can 
be implemented routinely. In this regard, in several recent studies 56,54 we have conducted 
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detailed investigations pertaining to both these approaches. The general conclusion drawn 
from these studies is that in cases where the AMC has proven useful, other approaches based 
on assumed probability distributions are also effective. In the cases considered in Ref. , 56 it 
is shown that in simple flows where the AMC can be employed, the class of PDF’s based on 
Johnson Edgeworth Translation (JET ) 57,58 can also be used. In fact, for the simple problem 
of binary mixing in isotropic turbulence - a standard test case - the solution generated by 
AMC 51,59 can be viewed as a member of the JET family. Furthermore, due to established 
similarities of JET with the simpler distributions belonging to the Pearson Family (PF) of 
PDF’s , 60 it can be argued that the PF can also be considered as a viable alternative. 

In turbulent combustion there has been widespread use of PF assumed PDF’s (e.g. Refs .; 61-66 
for recent reviews see Refs. 18,1,67 ). In most applications to date, this family has been used in 
the form of the Beta density of the first kind (Pearson Type I and Type II distributions). This 
is due to the flexibility of this density in portraying bimodal distributions. The capabilities 
of this density for the purpose of RANS have been studied in Refs . 52,56 ’ 54 According to these 
studies there axe some similarities between the PF and the AMC, as well as some distinct 
differences. As indicated before, both these methods are utilized in the context of a single- 
point closure. Therefore, in both cases the magnitudes of the moments up to second order 
must be provided externally. Also, both methods suffer from an inability to account for the 
shrinking bounds of the scalar domain as mixing proceeds. This is manifested by the failure 
of both closures in producing a correct evolution for the conditioned statistics of the scalar; 
namely, the conditioned expected dissipation and the conditional expected diffusion. This 
behavior could be troublesome and may produce significant errors especially in predicting 
the rate of reactant conversion in non-equilibrium flames. In their application for modeling of 
reacting homogeneous flows, both closures axe satisfactory for equilibrium flames regardless 
of the equivalence ratio . 52,54 However, the use of AMC is very difficult, if not impossible, for 
modeling of non- homogeneous flows regardless of the chemistry model. In such systems the 
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mapping transfer function must be evaluated numerically and in the non-equilibrium case 
the multivariate form of the PDF must be considered. These require further investigations 
before they can be implemented routinely. 51 In these cases the application of the PF is much 
more straightforward but obviously cannot be justified rigorously for all applications. The 
corresponding multivariate form of the Beta density is the Dirichlet distribution, 68-71 and its 
mathematical implementation is straightforward. 

Here, we report the results of our investigations pertaining the use of the Pearson family 
of PDF’s in LES of turbulent reacting flows. This investigation is based on a priori and 
a posteriori analyses of simple homogeneous and shear flows under reacting, non-premixed 
conditions. Both equilibrium and non-equilibrium chemical reactions are considered. A 
priori analyses are performed of equilibrium flames in both flow configurations. A posteriori 
analyses are performed only for the shear flow with non-equilibrium chemistry. In all cases, 
the LES generated results are appraised by comparison with data provided via DNS. 


2 Flow Configurations 

Two flow configurations are considered: (1) two- and three-dimensional homogeneous box 
flows, and (2) a two-dimensional spatially developing planar mixing layer. In both config- 
urations, a constant rate, non heat releasing chemical reaction of the type A + B — *• P is 
considered. 

The procedure in homogeneous DNS is similar to that of previous simulations of this type (see 
for example Ref. 1 ). The computational package is based on a modification of the computer 
code developed in Refs. 72-74 The code makes use of the spectral collocation algorithm with 
Fourier expansion functions 75-77 to approximate the spatial derivatives at N grid points in 
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each direction. It is now widely recognized that the collocation method is less expensive than 
other spectral variants ( e.g . Galerkin and Tau methods) in evaluating nonlinear terms. The 
aliasing errors inherent in collocation simulations can be effectively reduced by de-populating 
the large wave number interval (this can be accomplished by an isotropic truncation). To 
describe the temporal evolution, a time-splitting method is adopted in order to account 
properly for acoustic effects. The advancing in time is achieved using a third-order Runge 
Kutta scheme with a variable step determined from the Courant-Friedrich-Levy (CFL) sta- 
bility condition. The hydrodynamic field is assumed isotropic and is initialized in a manner 
similar to that in Refs. 73,78 In this setup, the energy contained in spherical shells is obtained 
applying an imposed energy spectrum function E(ic') = A; 4 exp( 2fc 3 /&o) (where repre- 
sents the integral scale of the flow) to the randomly distributed Fourier components of the 
scalar fields. The compressible energy content can be varied by means of a proper parame- 
ter related to the energy ratio between the two velocity modes (solenoidal and irrotational) 
resulting from the Helmholtz decomposition. No forcing mechanism at the low wave num- 
bers is applied; therefore, the turbulent field decays in time. Several other parameters are 
necessary to control the compressibility level in the turbulent field. These are the intensities 
of density and temperature fluctuations in addition to the root-mean-square (RMS) of the 
Mach number M t = The reacting scalar fields, assumed dynamically passive, 

are configured with smoothed square waves (error functions) in anti-phase, with the two 
species at the stoichiometric conditions. This is identical to the conditions considered in 
previous simulations of this type. 79,48,80 

The mixing layer configuration consists of two co-flowing streams traveling at different veloc- 
ities and merging at the trailing edge of a partition plate. 81-88 Two reactants, A and 5, are 
introduced into the high- and the low-speed streams, respectively. To expedite the formation 
of large scale vortices, low amplitude perturbations are imposed at the inflow plane of the 
layer. The flow field which develops in this setting is dominated by large scale coherent 
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structures. The numerical procedure is based on a hybrid finite-difference/pseudospectral 
scheme. A third order upwinding method is used in the streamwise direction and a spectral 
collocation method employing Fourier expansions is used for cross stream discretization 1 
Time discretization is done by the Adams-Bashforth technique. The computational domain 
is bounded by 0 < x < 326, —86 < y < 86, where 6 is the vorticity thickness at the in- 
flow. The highest resolution consists of 512 finite difference nodes and 256 collocation nodes. 
With this resolution, reliable DNS with a Reynolds number Re = 1000 and a Damkohler 
number Da = 10 (based on the velocity difference and the vorticity thickness at the inlet) 
are possible. 

In both flows, all the species (A, B , P ) are assumed thermodynamically identical and the 
fluid is assumed to be calorically perfect. The value of the molecular Schmidt number is set 
equal to unity in all simulations. It is also assumed that there is no trace of one of these 
species at the feed of the other one, i.e. complete initial segregation. With unity normalized 
mass fractions of the two species at their respective streams, a mixture fraction denoted by 
J can be defined: 18 




YAfajj) - yfl(x,,t) + 1 

n * 


Je [o,i]. 


(i) 


where Y{ denotes the mass fraction of species t. With total initial segregation the magnitude 
of the mixture fraction is unity in the stream containing A, and zero in the stream containing 
B. In non-reacting flows, the trainsport of this mixture fraction portrays the mixing process. 
In reacting flows, the variable J denotes a conserved “Shvab-Zeldovich” scalar variable. 18 
In the limit of an infinitely fast chemical reaction adl the statistics of the two reactants aire 


1 Within the past year, Geoffrey Shieh has modified this code by employing a sixth order compact param- 
eter finite difference scheme for discretization in the streamwise direction. We have obtained many exciting 
results pertaining to both DNS and LES. However, these results are not at a stage that can be included in 
this report. 
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directly related to those of the Shvab-Zeldovich variable: 32,89,18 ’ 47 

The assessment of the models is facilitated by direct comparison with DNS data. The 
resolution in DNS is dictated by the magnitudes of the physical parameters, with sufficient 
testing on the independency of the results to the number of grid points. All the subgrid 
statistics are constructed directly from the filtered DNS data bank. In a priori analyses the 
statistical properties of subgrid data are considered. In a posteriori analyses, the results 
predicted by LES on coarse grids are compared with those of filtered DNS. These analyses, 
correspond effectively to simulations with a homogeneous box filter in which the values of 
the filtered means are constant within the box. The ratio between the coarse grid and the 
fine grid is, therefore, a measure of the filter width. That is, it determines the scale at which 
the subgrid fluctuations are not accounted for deterministically. Consequently, the statistical 
behavior of the variable within the subgrid is dependent on the magnitude of the sample 
size there. If this size is too small, statistical analyses are meaningless. If it is too large, the 
essential features of large scale transport are masked. In a priori analysis of the box flow, a 
coarse grid consisting of 16 3 grid points are considered. In the shear flow, LES is performed 
on the domain with 64 finite difference grid points and 32 Fourier collocation points. 


3 Problem Formulation. I: DNS Framework 

Under the assumptions of ideal gaw and negligible body forces the set of nondimensional gov- 
erning equations (continuity, momentum and pressure) in Cartesian tensor notation, which 
relate the thermodynamical variables p (mass density), tij (velocity components), p (pres- 
sure), T (temperature) are given by: 
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( 2 ) 


dp dpuj 
dt dxj 


= 0 , 


dpuj djpujUj) _ dp 1 dajj 

dt dxj dii Re dxj ’ 


(3) 


?p + u + 7P ^i = 1 "U— * 

dt 3 dxj 1P d Xj MIPrRe dx } \ dxj Re 


(4) 


In these equations, the indices i and j are integer values adopting values 1, . . . , where 
£ = 2 in two-dimensional (2D) and £ = 3 in 3D simulations. The viscous stress tensor < 7 tJ 
and the viscous dissipation $ are defined, in order, as: 


°\j = P* 


( duj | duj 2 . dun 

dxi ) 3^ dxk u ’ 


(5) 


$ = /i* 




( 6 ) 


For a binary chemical reaction of the type A + B — *■ Products the above system is supple- 
mented by the equations of conservation of species: 


dpYg d(pY a Uj) 
dt dxj 


_L d_ 

Sc Re dij 



+ 


a = A,B 


(7) 


Here Y a denotes the mass fraction of species a, T is the binary Fickian diffusion coeffi- 
cient and u> a = —DapYjPpB indicates the chemical source term. The Damkohler number 
Da = poLokj/Uo, where kj is the rate constant of reaction and gives a measure of the speed 
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of reaction. In the present formulation, for simplicity, the diffusion coefficients are considered 
constant and the Soret-Dufour and pressure gradient diffusion effects are neglected. More- 
over, it is assumed that the Lewis number is unity and the chemical rate is independent of 
temperature. The system is closed with the equation of state for perfect gases. 


7 MlcP = pT. ( 8 ) 

The equations are scaled by reference values which appear in the expressions of the dynamical 
parameters: the characteristic Reynolds number Re = P 0 U 0 L 0 / po, the Prandtl number 
Pr = Cyfio/ko, the Schmidt number Sc = p 0 /T 0 and the reference Mach number. 


4 Problem Formulation. II: LES Framework 


In most cases of practical interest the motion on the order of the dissipation scale cannot be 
evaluated explicitly since the available computational resolution fails to meet the severe mesh 
requirements. In the homogeneous flow simulation, to insure a resolved flow simulation one 
has to consider a number of points in the range of iV 3 ~ 0.1 Re\ , while Re\ is the Reynolds 
number based on the Taylor microscale and RMS of velocity fluctuations 2 ). This implies 
that a 96 3 DNS can accommodate at most a modest Rex » 35. In order to overcome this 
limitation the governing equations have to be “altered” in such a way that the activity at 
the level of the unresolved scales is mimicked by a proper model, and only the large-scale 
fluctuations are explicitly computed. The rationale for this idea remains the fact that the 
primary momentum transport and turbulent diffusion is sustained by large-scale energy- 
containing eddies. Certainly a smoothing (“low pass” filter) operator achieves this goal, 
decomposing a given field into a resolved component and a residual component (also called 
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sub-grid field). 


It is customary to define for a flow variable /(x, i) the resolvable-scale component J through 
a convolution of / with a filter function 90 (7(x, t, A/): 


/(x, t) = [ G(x - x', t, A /)/(x', tJdV (9) 

where x = (z,-, i = and A/ is the filter width representing the smallest scale of the 

resolved field (it is desirable that A / corresponds to a wavenumber in the inertial range) and 
D is the domain of interest. Several types of filters are employed in previous attempts in 
LES: The Gaussian filter, the sharp-cutoff filter and the box filter. 90 The convolution with 
these filters is equivalent with a moving volume average and generates moderate subgrid 
fluctuations; therefore the modeling of the subgrid correlations is less critical than in the 
case of a fixed average. However, additional terms have to be modeled and an ill-posed 
Fredholm integral equation has to be solved in the necessary process of de-filtering. 

For the ease of implementation we employ the fixed average approach (used in a different 
manner in Ref. 20 by defining the “filter” function as: 


G(x,t,A f ) 


1 

0 , 


if 


1 - 


n sa 
otherwise 


A/ < X, < ttz! *A /; 

n sa 


( 10 ) 


with n = 1, Nsg represents the index of coarse grid points in each direction and tisg is 
the number of points per subgrid cell (box) side. Evidently the “filtered” variable is now 
defined only at discrete points. In the study of compressible turbulent flows it is convenient 
to prescribe a Favre filtered variable: 
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( 11 ) 



along with the corresponding decomposition of the instantaneous field 


/ = / + /", ( 12 ) 

which defines the SGS fluctuations /". The direct application of the filtering operator to the 
conservation of mass equation yields: 


&P . =0 

dt dxj 


( 13 ) 


where the commutativity of linear operators is exerted with the caveat that the density and 
the pressure can be meaningfully related only to a non-weighted filtered quantity. Similarly, 
the Navier-Stokes equations become 


dpui d(pUjUj) _ dp 1 daji d(-pUjUi) ■ - ... 

dt dij dii Re dxj dxj 

with 7 = pT. The evaluation of the nonlinear term pu t Uj gives: 


( 14 ) 


■pupTj = p(uiUj + u'/uj + Uiu'j + u"u"). 


( 15 ) 


With this relation the stress tensor is composed of a so-called “cross stress” Cij = —p(u"uj + 
l fpu'f) and the SGS Reynolds stress = -pu'-u". From this definition, it is apparent that in 
the present filtering procedure the cross stress is zero. Here, we employ the linear gradient 
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ii T!T ir 


class of isotropic models, which relate the SGS Reynolds stress to the large-eddy rate of 
strain Sij = ( dui/dij + diTj/dii ): 


Tii = -7>W< - = M(2 5,'j - |s«*y) - |jJ S6 ti (16) 


where the SGS Reynolds stress is separated into anisotropic and isotropic parts. The model 
developed by Speziale and co-workers is a compressible generalization of the Smagorinsky 
model which postulates that the eddy viscosity v t = C/jA^Iiy 2 , i.e. i/ t is related to the 
second invariant of the large scale strain tensor II 5 = Ski Ski- Although the isotropic part 
is relatively small compared to the thermodynamic pressure we prefer not to neglect it and 
use instead a model elaborated by Yoshizawa giving 2 = C/Aylfs. Cr and Ci are empirical 
constants to be determined. An alternative procedure is analogous to the “one-equation 
closure” in RANS. With this closure, the viscosity is expressed as i/< = C v Afs/% in terms of 
the SGS kinetic energy k = In this case 2 = k and Ci — 1.0. The term u : ) is 

treated like The filtered pressure equation becomes: 


dp _ dp 

dt + Uj dTj 


+ 7 P 


d£j 

dij 


M^PrRe dx 




+ 


1 d{-pu"T") 

7A/2, dij 


+ 


\i — 1 du*- dv> tf 


(H) 


[ 


The SGS velocity- temperature correlation Aj = — pu"T" is modeled based on the gradient 
transport approximation: 


A j = P 


vt df 
Pr t dij 


(18) 
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With this, the turbulent Prandtl number Pr t is taken as an empirical constant, implying 
the assumption that the transport of momentum and scalar quantities is characterized by a 
single length scale. Also, since the term u" can be interpreted as a SGS mass flux, we can 
write: 


~n Aj _ ^ dp 

j P Pdx / 

In the Smagorinsky model the eddy viscosity, vj, is related to the large scale strain rate. A 
suggestion made in Ref . 91 and discussed in Ref . 92 is to model the eddy viscosity in terms of 
the trace of the resolved vorticity field. Here, since the SGS moments up to second order axe 
required for the species field, an extension is made to utilize a one-equation hydrodynamics 
closure. The conjecture utilized for the molecular stress is also applied to the moleculax 
heat flux and dissipation. A modeling of the pressure dilatation term would require an 
additional transport equation of a subgrid quantity, for example the SGS density variance . 93 
Therefore, we neglect this term for simplicity. This is justified only when consistent decay 
occurs, that is after the adjustment between the thermal energy and TKE modes from the 
initial conditions is achieved. To complete the hydrodynamic system, an equation for the 
SGS turbulent kinetic energy is required: 



dpk d(pkuj) 
dt + dxj 


_d_ 

dx 


- [-(?<* + w - A«)] - ?«|| + 


^ dx, Re ji d Xj ia + 


1 —da,i 

-U; — 


dii Re * dxj ' 


( 20 ) 


The classical closure used for the SGS transport term Tj = —{pu'-k -f p'u'j — reads: 
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( 21 ) 


r,= 


(m J_\ dk_ 

V ft Re) dij' 


where ft is an empirical constant. The SGS TKE dissipation, c = contains a 

solenoidal part and a dilatational part, the latter being associated with the SGS viscous 
losses. A possible closure for the SGS dissipation including the compressible effects is: 


e = C ( [l+c d (M t SGS n (22) 

A/ 

Here C t is the SGS dissipation constant, c^ — 1.0 and Mf GS denotes the RMS of the SGS 
turbulent Mach number. The Favre filtering applied to the conservation of species equations 
demonstrates the coupling effects of ail the aerothermochemical processes in the density- 
weighted variables. With this, we have: 


&pY a d(pY a Uj) 
dt dxj 


_! a rrw^\ 

ScRedij \ dxj J dij 


a = A, B. 


(23) 


The Favre flux of species fluctuations :a represents the interactions between the two prin- 
cipal phenomena: the turbulence and the chemical reaction. Based on the previous assump- 
tions regarding the scalar fields, we have 






(24) 


where Set denotes the turbulent Schmidt number, also an empirical constant. The chemical 
source term establishes the rate of reactants consumption. This contains both large scale 
interactions and chemical SGS activity: 
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Z£ = -pDa(Y A Y B + YZYg). 


(25) 


The SGS covariance Y^Yg is of the same order of magnitude as the mean reaction term 
YaYb, consequently introduces an additional transport equation: 


gz n + mi Y B 5 >) 


i 


dt 


dxj 


ScRe dx 




C7WWE\ 

: r^f) 

„ WSSYg 


a(-KTO) 


ScRe dxj dxj 


+ u A Yg + u B Yjt + 


ScRe 


( r ‘ y » 


3^4 


5zj dxj 


A /W^ _ 

dxj { A dx, j dxj dx, 


(26) 


In the above equation, the last term on RHS is relatively small for flows with large Peclet 
numbers or if Y” is negligible. The gradient of the SGS correlation flux is transformed in a 
diffusion term expresses as: 


vS V u _-*< 


-KW =7 


S^ dxj 


(27) 


The simplest model for the SGS covariance dissipation c A b = * n terms of 

Y%Yg and the time scale of the turbulence r (the Smagorinsky model states that r = l/-y/TTs, 
whereas the one-equation model prescribes r = A//\/]b) is given by 


(-ab = C# 


*7*1 


(28) 


The remaining unclosed terms from this respectable set of equations axe adso among the most 
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difficult to model. Obviously, triple correlations like 0/4 Vjgf convey information from the high 
order subgrid statistics to the large scale field. These terms could be directly determined 
from the SGS joint PDF (if known) as the expectation of an arbitrary function M of the 
random variable ip = (Y A , Vb): 


M(rP) = f M(rP)P AB (iP)dxP . (29) 

Here Pab is the SGS Favre PDF obtained from the joint SGS PDF of density and species 

Pab{iP,p) 

Pab{4>) = = / P p AB{tp,p)dp. (30) 

Since the PDF of a random variable embodies the entire stochastic information pertaining 
to that variable, it is very tempting to consider a multivariate PDF of ail the dynamical 
variables for a complete description. A general transport equation can be derived in order to 
depict the evolution time. Unfortunately, the complexity of the full PDF method makes it 
very cumbersome for practiced purposes. 36 Restricting this approach only for the treatment 
of the chemical species still poses serious closure problems. 47,23 A less rigorous technique, 
but extremely facile, is to assume apriori the shape of the joint PDF which will represent 
the statistical behavior of the reactants within the subgrid. 

As indicated before, from the family of many different densities, the Pearson system of 
frequency curves and its multivariate generalizations have been found to display certain 
advantages, e.g. flexibility, over other possible candidates. In fact, the Beta density of 
the first kind yields a reasonable prediction for the probabilistic behavior of the conserved 
scalar variable in equilibrium flows (see the references cite in the Introduction). In general the 
Dirichlet distribution of the n species, maintaining the Favre convention, is given by; 70 . 6 ®.® 4 - 6 ® 
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n 

with the constraints 0j > 0, j = l,n; 13 0j < 1; >7* > 0, fc = l,n + 1. The n + 1 parameters 

i = i 

rjk can be evaluated if the values of an equal number of moments are supplied, therefore 
knowledge of at least one second order moment is necessary, a recognizable feature of single- 
point closures. It is straightforward to show that the product moments about the origin of 

n 

order r = $3 r j 3116 
j = i 
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where .S'* is the support of P(0) in the composition space. With this PDF, the mean of the 
mass fraction of species * is given by: 


0 ,' 


Vi 



k=l 


(33) 


Furthermore, the r,- + r^th order correlation (product moment about the mean) between 
species 0,- and 0j can be calculated as 


20 




/=! k - 1 



j ) (&)\'Fj) k — + ^ tlLilz IMgj ±Hin (^j + n — * — i) 

( E »?m + 1) • * • ( E Vrn + r, - / + r, - k - 1) 

m=l m=l 

(34) 


The variance of species * and the covariance of species i and j are, in order: 
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(35) 


In the particular case of our study n = 2, hence three parameters must be determined 
from known SGS moments. Our detailed investigations, some of which reported in Refs. 95,96 
indicate that there is no optimum choice between the moments up to the second order in 
the above setting, but two equally valid combinations 2 . These are provided by the subgrid 
averaged scalars Ya, Yb along with either the SGS scalar energy E# = Y'p + Yg 2 or the SGS 
covariance Y^Yg. In general we have: 


T}1 = E4A; T)2 = Yb A; j? 3 = (1 - Y A - Y s ) A; A = rji + t / 2 + rj 3 . (36) 

Also the unclosed terms in SGS covariance transport equation (which appear in the same 
form in the SGS scalar energy transport equation) are closed by means of the assumed PDF 
method in terms of known quantities as: 


u a Y: = Y a Y b { 1 - 2Y a ) 


A 

(A + 1)(A + 2) 


(37) 


2 Thi8 issue is being seriously investigated by George Sabini in a detailed study in which all the results 
generated by various assumed PDF’s are being compared with laboratory data. 
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The parameter A has the form: 


x _ Y a {1-Ya) + Yb{1-Yb) 

Yf + Yp 

when the scalar energy is used, and is given by 

A = + 1) (39) 

1 A X B 

5 Results. I. Homogeneous Box Flow 

5.1 DNS and a priori LES Analyses 

One objective of our direct simulations is to provide an extended data base for the pur- 
pose of SGS model validation. Although from a priori analyses one may extract a rather 
“pessimistic” conclusion about the performance of the models, it is the preferred manner to 
determine the best estimates for the unknown parameters appearing in the closures. In addi- 
tion, the ultimate test for a LES formulation is the comparison between a parallel evolutions 
DNS-LES. The most disturbing drawback of the single-point closures is the flow dependence, 
which can be alleviated to some extent by testing and fine-tuning under various conditions. 
The majority of the high-compressible turbulent flows include noticeable zones where the 
character of the flow shifts toward the quasi-incompressible regime. Hence it is required to 
endow the models with some flexibility. Finally, the assumptions made with respect to the 
magnitudes of neglected terms can be assessed by examining the balance of the equations. 

The present simulations, based on 256 2 (in 2D) and 96 3 (in 3D) collocation grids, are con- 
ducted to attain marginal levels of compressibility, that is either a weakly compressible 
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turbulence or towards a relatively high dilatational flow. As remarked before the initial 
conditions play a primary role in placing the simulation in the desired realm. ® The results 
of recent investigations at ICASE, NASA LaRC reveal that for 0 < M t < 0.3 and initially 
Prms * Trms < M t a quasi-incompressible turbulence statistics prevail, whereas for values 
of order the compressible TKE reaches the order of magnitude of the incompressible 

TKE. The set of initial parameters establishing the low compressible regime was M t = 0.2, 
Prms = 0.01, T RMS = 0.01, x = 0.01, where x = £?!(£? + S 1 ). To produce stronger 
compressibility effects, computations are performed with the initial conditions M t = 0.6, 
Prms = 0.1, Trms = 0.1, x = 0.2. The box Reynolds number is taken in the range 100 - 250 
(the lower value corresponds to the high compressible case). 

The above mentioned effects can be evidenced qualitatively through visualization of relevant 
quantities and reflected quantitatively in the time history of various ensemble averages. An 
illustrative contour plot of the thermodynamic pressure, presented in Fig. 1 allows the identi- 
fication of small compression waves (shocklets) within the high compressible turbulent field. 
For extensive phenomenological descriptions we refer the reader to our previous works.®® 
From the balance of the TKE equation (Fig. 2) for the 3D compressible case (with initial 
R&\ — 12) it can be pointed out that the pressure dilatation term creates strong interference 
in the early phases of the energetic decay (all terms are obtained by Favre ensemble average). 
Moreover, the dilatational dissipation represents a considerable part of the total dissipation 
(about the same percent of compressible energy injected into the flow). These contributions 
have to be taken into account by means of external models. 

The evolution of the reacting scalars is a reflection of the development of the flow. Figure 3 
presents the contour plot of species A in a 2D compressible simulation with Da = 10. The 
figure shows that the initially segregated structure has been contorted by turbulence and 
the species are depleted by chemical reaction. The influence of the source term on the decay 
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of the second order moments can be appraised by Fig. 4 where we present the balance of 
the transport equations for the reactants covariance and the variance of species A. It can be 
noticed that the contribution of the reaction term is significant, affecting the rate of decay 
set up by the scalar dissipation. Figure 4(c) justifies the neglect of the last term in Eq. (26) 
based on small values of Y%. 

In order to evaluate the performance of the Dirichlet joint PDF an a priori analysis is 
conducted on the data obtained from high compressible 3D non-equilibrium chemistry case 
with Da = 5. A striking feature displayed in the statistics of the reacting scalars is that 
the density corrupts the time history of all the moments except the second order ones which 
exhibit essentially the same profiles under both ensemble and Favre ensemble averaging 
procedures (Figs. 5(a)-(d)). For practical purposes, a reacting flow is completely determined 
if the scalar’s moments up to the second order are known, hence the particular case of 
binary systems necessitates five quantities. Although the Dirichlet density is not able to 
accommodate five parameters, it is preferred for its flexibility and tractability. Obviously, 
an optimum choice has to be made based on the best approximation of the scalar covariance 
and of third order moments needed to be closed in the transport equations of the lower order 
moments. The examination of the Figs. 6(b) and 7(b) gives more credit to the means-scalar 
energy parameterization, but a careful reconsideration of the fact that the source term in the 
species conservation equation is under the direct influence of the scalar covariance evolution, 
shifts the option towards the means-covariance triad. Figure 6(a) enforces whereas Fig. 7(a) 
weakens this choice. 

More insight can be acquired exploring the evolving shape of the scalar joint PDF (based 
on the fine grid statistics) as resulting from DNS. Figures 8, 9 and 10 portray this evolution 
in the composition space. At initial times when the scalars A and B are still isolated (to be 
rigorous, there are overlapping regions due to the initial conditions-error functions) the PDF 
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is roughly composed of two impulse functions at the higher bound. Later, the PDF breaks 
into a wave-like surface, which travels to the lower bound. Towards the completion of the 
reaction the PDF is converted to an unique impulse shape situated at tp = 0. An examination 
(Fig. 11) of the subgrid joint PDF’s exposes skewed shapes and similar trends as described 
before (the resolution is insufficient to obtain smooth surfaces). From a purely geometric 
standpoint, the PDF can be depicted as a family of Beta functions (in vertical planes) having 
as support a hyperbola in the horizontal plane (Fig. 12) which is the contour representation 
of Fig. 9. Unfortunately, the Dirichlet density does not possess these geometrical features as 
portrayed in Fig. 13 (this case corresponds to the parameterization with the scalar energy E+ - 
the one based on covariance produces a similar shape). This mismatch is evidently translated 
into statistical disagreements. In an effort to overcome this deficiencies we propose a PDF 
(Fig. 14) which holds sufficient degrees of freedom and also behaves similarly to the joint 

PDF of reacting scalars: 


^ w + v-r'w. + - * - * - ( 4 °) 

where C N is the normalization constant. For the case considered here, this PDF is very 
accurate. However, there are two severe drawbacks which does not allow its use for general 
applications. First, the extension of this PDF for multiple scalars is not trivial. Second, it is 
not possible to relate the parameters of this PDF to the moments (or cross moments) of the 
random variable. While the first problem can be potentially overcome with availability of 
DNS data, the second problem is the compelling reason for abandoning this PDF for LES. 

The data generated from DNS is now employed to validate the subgrid models and to deter- 
mine the unknown model constant. In both 2D and 3D, a coarse grid is superimposed on the 
fine grid (, e.g . 256 2 - > 32 2 and 96 3 - > 16 3 ) and all the relevant quantities are determined: 
averages, fluctuations and correlations (stresses and fluxes). The second order information 
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which is not available in LES is compared against the proposed closures. The inspection and 
statistical analysis of the scatter plots of the exact quantities vs. the predicted values provide 
qualitative and quantitative assessments. The correlation coefficient between the exact SGS 
stress and its closure is defined as: 


< (E- < E >){M— < M >) > 

Cem ~ <{e-<e >y >V2< (m- < m >y >*/2 


(41) 


where the brackets denote the ensemble average. This correlation coefficient provides a 
proper mean of diagnosing the reliability of the closures. Examining the time histories of the 
correlation coefficients extracted for each SGS variable from 2D and 3D simulations, it can be 
concluded that the Smagorinsky closure and the one-equation model attain similar degrees 
of reliability even in compressible regimes. This point is illustrated in Fig. 15(a)-(d) for the 
xy SGS stresses, SGS species and covariance fluxes in x-direction and scalar dissipation. The 
degrees of correlation are relatively poor except for dissipation closures. A large variety of 
algorithms can be employed in the calculation of the model constants at various levels. We 
made determinations only at the scalar level using three techniques: mean matching, RMS 
matching and linear regression. The matching method involves the computation at each time 
step the ratios between the ensemble averages and between the RMS values, respectively: 


„ _ <E> „ _<(£-< E >y > 1 / 2 

A < M > ' ( " RMS ~ < (M- < M >) 2 > x / 2 ’ 

The linear regression method is applied to the scatter data obtained from time steps within 
a large eddy turnover time. The mean matching method produces variations of the ratio 
over several order of magnitudes. This reflects the deficiency of the models (Fig. 16). A 
smoother evolution is obtained through RMS matching (Fig. 17). The decisive estimation is 
provided by the linear regression based on a larger sample of data (Fig. 18) which displays 
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indeed relatively strong dispersion from a linear variation. Through this final procedure we 
fixed the values of the constants as follows: For the Smagorinsky model 


C R = 0.16; C/ = 0.008; C+ = 5.9; 


(43) 


and for the one-equation model 


C u = 0.034; C t = 0.78; C, = 2.1. (44) 

From the scatter of the species flux a turbulent Schmidt number was estimated at Sc< = 0.77. 
Also, this value was assigned to the turbulent Prandtl number. 


5.2 LES: A posteriori Analysis 

The effective performance of the SGS models have to be proved by conducting actual LES 
with the same dynamical parameters as the peer DNS. In order to start the LES from an 
equivalent initial state the DNS fields are filtered and introduced into the coarse grid. The 
comparison between parallel simulations is achieved at each time step by filtering the DNS 
output. 

The first test of the LES framework presented here is related to the fact that the assumed 
PDF has a considerable influence over the reacting scalar evolution. We used the results 
of 3D DNS of a binary mixing (Da = 0) field in a weakly compressible turbulence as the 
prime test for LES. In this case the transport of the SGS covariance decoupled from the 
other equations. The obliteration of the source terms in the chemical equations also removes 
the weight of the assumed PDF. In Fig. 19 the mixing structures are presented by means 
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of contour plots of the species field. These results tire presented at one half turnover time. 
There is an overall agreement between the two pictures, the large scale structures being 
predicted fairly accurately. However, the small scale structures appear to be different in 
some areas. Comparisons can be made in a global sense considering various moments of 
the LES scalars based on ensemble average. The resulting correlation (connected to the 
concept of unmixedness), the variance, the kurtosis and one of the cross third moment show 
satisfactory concordance with the DNS results (Fig. 20). A more stringent test stemming 
from the contrast at the subgrid level divulge inadequate evolutions of the second order SGS 
moments. Figures 21 and 22 depict the time evolution of the SGS variables at different 
locations (starting with different initial values). Although the SGS mean species evolve 
in group, there are visible discrepancies in the predictions of the covariance and turbulent 
kinetic energy, as a consequence of the inferior quality of the models. 

The non-equilibrium chemistry case was investigated for a relatively high rate of reaction 
Da = 5. In the present situation, a correct appreciation of the output from LES must be 
weighted with respect to the degree of success exhibited by the PDF closure. Namely, after 
examining Figs. 23 and 24, the results could be considered “satisfactory.” The predicted 
ensemble averages of the product and of the scalar A tend to give to the k — A model an 
advantage over the Smagorinsky closure. Further examination of the higher order moments 
enforces this opinion. 


6 Results. II. Spatially Developing Planar Shear Flow 

In employing assumed PDF methods in an actual LES, the filtered mean values and SGS 
moments of the scalar variable must be provided without resorting to DNS. Here, LES re- 
sults obtained by the full Smagorinsky /PDF closure are presented to examine the overall 
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performance of the hybrid model. This examination is conducted for the mixing layer since 
it provides a more stringent test for model assessments. The formulation presented above is 
in a closed form with the specification of model constants appearing in the LES transport 
equations. Based on our investigations, it was clear that the magnitudes of these constants 
are different from those in the homogeneous flow simulations. Therefore, a nominal para- 
metric study was done to re-estimate the values of these constants for the mixing layer 
simulations. From this study, the magnitudes of some of these constants were found also 
different from those typically used in RANS. 97 The results presented in the next section are 
based on C + = 0.15, Set = Cjt = 1.0, C v — 0.01, and C t = 0.5. No attempt was made to 
determine the optimized values of these constants. 

In Fig. 25, cross-stream profiles of the instantaneous streamwise velocity component are 
shown at several downstream locations. This figure shows that the large scale velocity field 
is predicted well in LES as the results compare reasonably well with those of filtered DNS. 
The comparison of the subgrid scale TKE as predicted by the model, with the TKE obtained 
from the filtered DNS are shown in Figs. 26-27. These figures also indicate a relatively good 
agreement between the LES predictions and the DNS filtered results. The exception is at 
far distances from the inlet, where it is shown that the TKE equation severely underpredicts 
the subgrid turbulence level. This trend has also been observed in the simulations of Ref. 22 
and may be due to deficiencies in modeling of the Reynolds stresses in terms of the strain 
tensor alone. Inclusion of other terms such as the rotation tensor and/or products of the 
strain rate and rotation tensors may be necessary to improve the predictions of TKE. 93,98 

Figure 28 portrays the LES and DNS contours of the filtered mass fraction of species A in the 
non-reacting case. The agreement is generally good and the two results axe relatively close. 
The difference is due to the effects of the gradient diffusion closure which results in smearing 
of the small scale features within the large scale coherent structures. The lack of agreement 
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becomes more pronounced for the covariance. This is shown in Fig. 29 where contour plots of 
the SGS species covariance are presented. Despite this difference, it is useful to examine the 
model’s performance in the reacting shear flow LES. In Fig. 30, contour plots of the filtered 
product species mass fraction obtained from LES and DNS are presented. Figure 31 shows 
an instantaneous product thickness distribution obtained from the filtered DNS and LES 
data, including also some LES results without the PDF model. These figures suggest that 
LES captures the large scale feature in accord with DNS. The figures also suggest that the 
inclusion of SGS fluctuations gives results which axe closer to DNS data than those predicted 
by the mean chemistry model. However, the PDF method still needs to be improved further 
in order to exhibit a better predictive capability. 

The primary reasons for the discrepancies above are attributed to errors associated with: 
(1) the estimation of the covariance, and (2) the shape of the PDF. We feel that with the 
simplified reaction mechanism considered in this flow , the first factor is more important. This 
is indicated in Fig. 32, where contour plots of SGS covariance are presented. Again, while the 
LES and filtered DNS results show similar trends, the agreement is not very good (see Fig. 33 
for a quantitative comparison of the covariance thickness). These two figures, along with Fig. 
29 highlight the modeling deficiencies associated with the covariance transport equation for 
LES applications. Therefore, while this equation has been met with some success in RANS, 
its use here does not provide acceptable predictions. The same is true for the SGS kinetic 
energy as shown in Figs. 26-27. In this regard, it must be noted that these problems do not 
vanish in PDF approaches based on a transport equation for the single-point PDF (AMC, 
C/D, LMSE, etc). That is, the first two SGS moments must be determined by a reliable 
subgrid closure. The problem may be somewhat alleviated by considering the transport 
equations for SGS correlations." This is very challenging and computationally demanding. 
Some saving in computations can be made by approximating the transport terms in the 
scalar flux equation in a way analogous to what is done in algebraic Reynolds stress closures 
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(see Refs. 100-102,93 ). The performance of these more elaborate models are presently being 
assessed for the modeling of SGS covariance and TKE. 


7 Some Remarks 

During the past century, dating back to the early pioneering contributions in statistics and 
biometrics , 60,58 ’ 57,103 the construction of “appropriate” PDF shapes has been a subject of 
broad interest within the statistics and (old) biometrics research communities. The out- 
growth of these contributions has been useful to investigators in other branches of physical 
science. In fact, in classical turbulence research, statistical methods 104 of one form or an- 
other had been the primary means of dealing with turbulence and its “random” causes and 
effects . 105-107 In statistical modeling of turbulent reacting flows, the use of PDF methods 
has proven particularly useful. This is especially true if the PDF (or joint PDF) of scalar 
quantities is considered. The preferred method is to obtain the PDF by means of its trans- 
port equation. This method is a subject of ongoing research 23 and there are still a number 
of questions in regard to its suitability for practical applications. LES can be viewed as an 
example of such an application, since it is viewed as more of a potential engineering tool 
than a robust scientific tool. With this view, assumed PDF methods are advocated as the 
present method of choice, at least in the context of the simple flows considered here. 

In this work, an attempt has been made to borrow knowledge from the statistics literature in 
order to “presume” an appropriate PDF which performs reasonably well in LES. In simple 
flows of research interest such as a homogeneous box flow or a parallel shear flow this is 
a viable approach in light of the richness of literature on their behavior (see Ref . 86 for a 
recent review). Consistent with the current state of affairs in RANS, these PDF’s are all 
considered in the context of a single point. This implies that all the moments up to the 
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second order (either Reynolds moments or SGS moments) have to be provided externally. 
Based on our earlier investigations along these lines, 62 ’ 54 ’ 5 ® the use of the Pearson family 
of PDF’s appears to be most practical. This suggests a Beta density of the first kind for 
equilibrium (and frozen) chemistry, and the Dirichlet density for non-equilibrium chemistry. 
In the specific cases considered, the use of the Beta density provides a reasonable model 
for SGS probability description of a conserved scalar variable. However, its actual use in 
LES requires the knowledge of the filtered mean and the SGS variance. This is currently 
troublesome in that conventional turbulence closures do not seem to work well for predicting 
second order SGS moments. The same problem exists in LES of non-equilibrium flows 
and manifests itself in the inability to accurately predict SGS scalar covariance (and/or SGS 
scalar energy). Moreover, it is not clear which of the sets of low order moments are to be used 
to parameterize the multivariate PDF. This drawback is, however, deemed less troublesome 
in that the Dirichlet density can be replaced with other members of the joint Beta family 
which are constructed with a higher degree of freedom. Again, the statistics and biometrics 
literature can be helpful for this purpose, and this issue is currently under investigation. 

At this point, it is useful to make some remarks in regard to some situations where assumed 
PDF methods are not very successful. These correspond to cases where an obvious choice 
for the PDF can not be made. For example, in highly non-equilibrium flames with complex 
multi-step multi-scalar chemistry with strong temperature dependent reactions rates. In 
LES of such flows, the application of the PDF transport equation may also be difficult. 
However, other less computational intensive closures may prove useful. For example, the 
“laminar diffusion flamelet model”, 108 ’ 641 109> 110,18 or the “conditional moment method ” 111-113 
may provide some alternatives. An appraisal study of the performance of these models in 
RANS of turbulent mixing layers has been reported in Ref . 114 The first model is applicable 
to cases where chemical reactions occur in a narrow region near the flame surface. For 
s uffi ciently fast chemistry, but with finite values of the Damkohler number, the flame is 
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located somewhere within the (unresolved) subgrid. The actual LES with this model requires 
the construction of a flamelet library within the subgrid by which the reacting scalar values 
can be approximated. The approximation requires the knowledge of the filtered mean values 
of the mixture fraction ( e.g. the Shvab-Zeldovich variable), and its SGS rate of dissipation. 
Therefore, modeled transport equations are needed to approximate these quantities. How 
effectively the dissipation field can be modeled and how accurate the whole procedure would 
be, needs to be determined. The second approach can be considered somewhere between 
the PDF approach and traditional moment methods. In its implementation for LES, the 
averages of the reacting scalar field are defined “conditionally” on the values of the mixture 
fraction. The argument in support of this closure is that such conditional statistics portray 
less scatter than their unconditional counterparts. Therefore, their treatment may be easier. 
However, an extra dimensionality (associated with the domain of the mixture fraction) is 
involved. Also, the actual application of the model requires the input of the SGS conditional 
expected dissipation. This could be very difficult since the behavior of this dissipation is less 
understood . 53 Furthermore, as discussed in Ref ., 114 the use of the model can be recommended 
only for determining the first order moments in flames near equilibrium. This is acceptable 
if the mean compositional structure in near-equilibrium flames is desired. Otherwise, the 
usage of the method is complicated in RANS and even more complex in LES. 


8 Further Information 

As indicated in the abstract of this report, here we have only discussed the results of 
our efforts pertaining to LES of turbulent reacting flows. Also, as indicated in the ab- 
stract, the scope of the research conducted under this grant is much more broad - perhaps 
much more than that can be described in a single report. No attempt is, therefore, made 
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here to describe all our contributions. However, for those who are interested we refer to 
Refs . 115,80,52,116,117 ’ 56,118 ’ 88 ’ 96 Copies of some of these papers are provided in Appendix I 
through Appendix VII. 

Mr. Craig J. Steinberger is the other Graduate Research Assistant supported partially by 
this Grant. For a summary of his contributions, please see the appendices. 
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Figure Captions 


Figure 1: Pressure contours for the high compressible case. 

Figure 2: Balance of turbulent kinetic energy transport equation. 

Figure 3: Mass fraction contours of reacting scalar A, Da = 10. 

Figure 4: Balance of transport equations for: (a) Reactants covariance, (b) reactant A 
variance, (c) time evolution of the mean magnitudes of the scalars and the scalar Favre 
fluctuations. 

Figure 5: Temporal variations of the scalar statistics under ensemble and Favre ensemble 
averaging: (a) Mean, (b) covariance, (c) variance, (d) third order joint moment. 

Figure 6: Temporal variations of the (a) second and (b) the third order Favre cross moments, 
using the means-scalar energy parameterization. 

Figure 7: Temporal variations of the (a) second and (b) the third order Favre cross moments, 
using a means-covariance parameterization. 

Figure 8: Snapshot in the composition space of the scalar joint PDF after 1/10 turnover 
time. 

Figure 9: Snapshot in the composition space of the scalar joint PDF after 1/2 turnover time. 

Figure 10: Snapshot in the composition space of the scalar joint PDF after one turnover 
time. 

Figure 11: Snapshot in the composition space of the scalar subgrid joint PDF after one 
turnover time. 

Figure 12: Contour representation in the composition space of the scalar joint PDF after 
1/2 turnover time. 

Figure 13: Dirichlet distribution with the means-scalar energy parameterization after 1/2 
turnover time. 

Figure 14: Proposed joint PDF after 1/2 turnover time. 

Figure 15: Temporal evolution of the correlation coefficients (between exact value and clo- 
sure) for SGS quantities: (a) xy stress, (b) species x-flux, (c) covariance flux, (d) scalar 
dissipation. 

Figure 16: Temporal evolution of the mean ratio (exact/closure) for SGS quantities: (a) xy 
stress, (b) species x-flux, (c) covariance flux, (d) scalar dissipation. 
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Figure 17: Temporal evolution of the RMS ratio (exact /closure) for SGS quantities: (a) xy- 
stress, (b) species x-flux, (c) covariance flux, (d) scalar dissipation. 

Figure 18: Scatter plots of the exact versus k — A closure for SGS quantities: (a) xy-stress, 
(b) species x-flux, (c) dissipation, (d) scalar dissipation. 

Figure 19: Mixing structures in a binary reacting field after one turnover time: (a) DNS, (b) 
LES. 

Figure 20: Temporal variation in the global statistics of the SGS scalars: (a) covariance, (b) 
variance, (c) kurtosis, (d) third order cross moment, Da = 0. 

Figure 21: Temporal variation of the SGS variables (location 1): (a) mean of scalar A, (b) 
mean of scalar B, (c) covariance, (d) turbulent kinetic energy, Da = 0. 

Figure 22: Temporal variation of the SGS variables (location 2): (a) mean of scalar A, (b) 
mean of scalar B, (c) covariance, (d) turbulent kinetic energy, Da = 0. 

Figure 23: Temporal variation in the global statistics of the SGS scalars: (a) mean product, 
(b) mean scalar A, (c) super-skewness, (d) fifth order cross moment, Da = 5. 

Figure 24: Temporal variation in the global statistics of the SGS scalars: (a) covariance, (b) 
variance, (c) kurtosis, (d) third order cross moment, Da = 5. 

Figure 25. Cross-stream variations of the of streamwise filtered velocity at several streamwise 
locations of the mixing layer, (a) x = 56, (b) x = 20 6, (c) x = 25 6. 

Figure 26. Plot of subgrid scale turbulent kinetic energy contours, (a) LES, (b) DNS. 

Figure 27. Cross-stream variations of the of subgrid kinetic energy at several streamwise 
locations of the mixing layer, (a) x = 56, (b) x = 206, (c) x = 256. 

Figure 28. Plot of subgrid averaged species A mass fraction contours for non-reanting mixing 
layer, (a) LES, (b) DNS. 

Figure 29. Plot of subgrid scale covariance contours for the non-reacting mixing layer, (a) 
LES, (b) DNS. 

Figure 30. Plot of product species mass fraction contours in the reacting mixing layer with 
Da = 10 (a) LES, (b) DNS. 

Figure 31. Streamwise variation of the instantaneous product mass fraictions thickness for 
the reanting mixing layer with Da = 10. 

Figure 32. Plots of subgrid scale covariance contours for the reacting mixing layer with 
Da = 10. (a) LES, (b) DNS. 
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Figure 33. Streamwise variation of the instantaneous covariance thickness for the reacting 
mixing layer with Da = 10 
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Abstract Results are presented of direct numerical simulations of a two-dimensional temporally dev elop- 
ing high speed mixing layer under the influence of a second-order non-equilibrium chemical reaction of the 
type A + B -> Products -j- Heat. Simulations are performed with different magnitudes of the convective 
Mach number and with different chemical kinetics parameters for the purpose of examining the isolated 
effects of the compressibility and the heat released by the chemical reaction on the structure of the layer. 
A full compressible code is developed and utilized, so that the coupling between mixing and chemical 
reaction is captured in a realistic manner. 

The results of numerical experiments indicate that at the initial stages of the layer’s grow th, the heat 
release results in a slight enhanced mixing, whereas at the intermediate and the final stages, it has a reverse 
influence. The effect of compressibility is the same in all stages of the development; increased compressibil- 
ity results in a suppressed mixing and, thus, in a reduced reaction conversion rate. Mixing augmentation 
by heat release is due to expansion of the layer caused by the exothermicity, and mixing abation is caused 
by suppression of the grow th of the instability modes due to increased heat release and/or compressibility. 

Calculations are performed with a constant rate kinetics model and an Arrhenius prototype, and the 
results are show n to be sensitive to the choice of the chemistry model. In the Arrhenius kinetics calculations, 
the increase of the temperature due to chemical reaction is substantially higher than that of the constant 
rate kinetics simulations. This results in a more pronounced response of the layer in all stages of the growth, 
i.e., an increased thickening of the layer at the initial phase of grow th, followed by subdued thickening at 
later stages. 

Key words: Supersonic combustion. Turbulent reacting flows. Direct numerical simulations, Compressible 
shear flows. 


1 INTRODUCTION 

Within the past five years, there has been a renewed interest in the analysis of high 
speed flows, both reacting and non-reacting. This interest has been, to a large extent, 
motivated by the need for the development of airbreathing vehicles capable of cruising 
at hypersonic speeds. This need in view of the goals of U.S. government in construct- 
ing such vehicles in the near future has resulted in a dramatic increase in research 
activities (both fundamental and applied) in the area of supersonic combustion. Of 
particular interest are the reactive systems of supersonic combustion ramjet (scramjet) 
engines (Waltrup and Billing, 1986; White el al., 1987; Drummond, 1991 ; Drummond 
et al ., 1991 Givi and Riley, 1991). According to recent surveys of airbreathing 
technology for hypersonic applications, in the combustor of a scramjet engine the fuel 
is injected at near sonic conditions into a supersonic air stream and combustion occurs 
in a mixing zone inside the combustor. The analysis of such a high temperature 
turbulent shear flow is complicated by the interaction between fluid mechanical and 
chemical effects, where the reactants first mix and then combine to release heat which 
in turn alters the structure of the flow. 

A basic understanding of the complex physical phenomena associated with high 
speed combustion requires a detailed knowledge of the effects of mixing and chemical 
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reactions, and the coupling amongst them. Analyses must he done first in a simple 
configuration so that the effects of the important physical phenomena can he some- 
what isolated and studied in detail. After the establishment of a basic understanding, 
it is then possible to make use of the knowledge gained for the purpose of design and 
for optimization of the actual applied processes. 

The exact processes of mixing and non-equilibrium chemical reactions in an actual 
scramjet engine arc far too complex to be investigated in a single study. Within the 
past decade, there has been much effort in research activities focusing on different 
aspects of high speed combustion, and the scientific community in this area has 
benelited substantially from the knowledge gained by such investigations (see Drum- 
mond, 1991 for a recent review). These activities have been concentrated on both 
analytical/computational and experimental efforts, and the results obtained to dale 
have been invaluable in providing a better understanding of the intricate details of the 
reacting compressible transport. The outcome of these analyses, together with the 
findings of ongoing and future scrutinies will undoubtedly result in an even better 
comprehension of these complicated phenomena. 

In these efforts, the results of which are reported here, we have concentrated on 
some specific issues regarding the transport of high speed chemically reacting flows. 
Namely, some effects of compressibility and heat release are the major issues under 
consideration. The approach is based on a computational effort, and among the tools 
currently in use, we have selected one based on direct numerical simulations (DNS). 
This approach, belonging to the class of “model free” methodology (in the terminol- 
ogy adopted by Givi, 1989), has proven quite useful for investigating reactive tur- 
bulent flows and its implementation for high speed combustion is deemed plausible. 
This has been convincingly evidenced by the outcome of several recent contributions 
in non-reacting high speed simulations (Lele, 1989; Mukunda et ai, 1991; Sandham 
and Reynolds, 1989; Soetrisno et al., 1988; Tang et ai, 1989). 

Since this is a focused (and somewhat a “narrow-minded") investigation of a very 
complicated phenomenon, we have selected a simple system in which each of the 
desired effects can be optimally elucidated. The system adopted is composed of an 
unsteady two-dimensional shear layer under the influence of a non-equilibrium 
irreversible chemical reaction. Due to limitations of the computational methodology, 
many simplifying assumptions regarding the physics of the (low and the transport of 
the aerothermochemical variables are made. These idealizations make the problem 
computationally more tractable, without disallowing the explication of desired physi- 
cal phenomena. 

In the next section, a description of the problem is provided together with a 
presentation of the non-dimensional grouping of the variables characterizing the 
physics of the problem. In Section 3, the results of simulations are given with a 
discussion on the actual physics portrayed by such results. The conclusions of this 
study are summarized in Section 4. 


2 DESCRIPTION OF THE PROBLEM 

The problem under study is a compressible two-dimensional temporally evolving 
mixing layer (Figure 1(a)). Within the “temporal" framework, the direction of the 
flow on the top stream is toward the right with a velocity U , . The stream on the 
bottom half of the layer flows to the left w'ith the same speed ( — U , ). The region of 
shear created by the velocity difference across the layer resembles that formed in a 
laboratory “spatially developing” layer (Dimotakis, 1989), and the temporal layer 
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FIGURE l (a) Schematic diagram of a temporally evolving mixing layer, (b) The grid. 


exhibits many features similar to those of a spatially developing one; namely, the 
formation of large scale coherent structures and their subsequent influence on the 
outcome of chemical reactions (Givi, 1989). The exception is the asymmetric mechan- 
ism of the mixing procedure in the spatial case (Dimotakis, 1986, Ghoniem and Ng, 
1987, Grinstein et al., 1986, Korczak and Hu, 1987, Lowery and Reynolds, 1986 and 
Givi and Jou, 1988 among others; see Givi, 1989 for a review) which cannot be 
captured in temporal simulations. Since this aspect of shear mixing is not the subject 
of the present investigation, provides a proper justification for the simulations in the 
format prescribed here. Further justifications of simulating temporally developing 
mixing layers, is provided in many analogous previous works (e.g., Claus, 1986, Givi 
et al., 1986, Metcalfe et al., 1987, Riley et al., 1986, McMurtry et al., 1989 among 
others). 

The speed in both streams is large enough so that compressibility effects cannot be 
ignored. A parameter which characterizes this effect is the magnitude of the free 
stream Mach numbers. For unity temperature and density ratios of tree stream flows, 
the Mach number is the same in two streams and is the same as the “convective Mach 
number, M c " identified and widely utilized in characterizations of spatially develop- 
ing compressible mixing layers (Papamoschou and Roshko, 1988, Dimotakis, 1989). 

The reacting species are introduced into the layer at the free streams. The chemical 
reaction occurring within the flow is idealized to a simple irreversible second-order 
form of A + B — ► Products + Heat. Reactant A is introduced on the top stream, 
and reactant B on the bottom stream. The two species are assumed to have identical 
thermodynamical properties, and the mass diffusion is assumed to obey the Fick's 
law. Moreover, the assumption of a calorically perfect gas is made with a constant 
ratio of the specific heats. 

The flow field is initialized with a hyperbolic tangent streamwise velocity with a 
specified initial vorticity thickness (<5 W | 0 ). At the initial time, species A is assumed to 
cover the top half of the layer, and species B the bottom half. The initial temperature, 
pressure and the density are assumed uniform throughout the field and the size of the 
computational box in the two directions (L v , L, ) is large enough to accommodate for 
the growth of the layer. In all the simulations LJ(SJ 0 ) — LJ(6 J 0 ) 5: 100. 
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The chemical reaction between l ho two species is characterized b\ the kinetics 
mechanism employed in the simulations. It is assumed that the rate ol reactant 
conversion is given by a single-step model ol the type: 


to — K t C A C H , ( 1 ) 

where C 4 and C H denote the concentrations of the two reacting species, and arc 
assumed equal at free streams, i.e.. C, , = C H< = C , . K, is the parameter measuring 
the extent of reaction. This parameter is normalized to form the definition of the 
Damkohler number: 


Da, 


*/Cl 
C, I<\ ; 1(1 


( 2 ) 


Two chemistry models are employed in the simulations, one with constant rate 
kinetics (i.e. constant K,). and the other with an Arrhenius model in which K , varies 
with the temperature, represented by the following equation: 


K, = 



( 3 ) 


where T represents the temperature normalized with respect to its free stream value, 
T f , Aj is the pre-exponential factor, and Ze is the Zeldovich number defined by: 


Ze 



(4) 


Here, E is the activation energy, and R is the universal gas constant. In the Arrhenius 
kinetics calculations, the pre-exponential factor ^ replaces K, in the definition of the 
Damkohler number (Eq. 2). 

Combustion exothermicity is measured by the energy liberated by the chemical 
reaction AH 0 . The magnitude of this energy is parameterized by a non-dimensional 
heat release parameter Ce defined by: 


Ce 


-AH 0 
C t ,T x ’ 


(5) 


where C., is the specific heat. In this form, Ce — 0 corresponds to a non-heat releasing 
chemical reaction. 

The governing equations of the compressible reacting flow field are solved by a new 
version of the SPARK computer code originally developed by Drummond (1987, 
1988). This new version, in comparison to the original form, has the capability of 
employing several high order finite difference schemes for solving the governing 
transport equations of compressible reacting flow fields. The numerical algorithm 
employed in this code is tailored for full compressible simulations, so that all the 
important aspects of compressibility (Chu and Kovasznay, 1958; Kovasznay, 1957) 
are captured. In its current form, the code is capable of implementing several ad- 
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vanced numerical algorithms lor this purpose (Carpenter, 19X8, 1989, I 1 ) 1 )!)). Among 
the available options, we selected to utilize the scheme suggested b> Goulieb and 
Turkel ( 1976). This scheme provides a fourth order accuracy in evaluating the spatial 
derivatives and a second order truncation in the temporal discret/ialion. The results 
of some rigorous tests by Carpenter (1990) and Steinberger (1991) show that with 
proper resolution, this algorithm is comparable with schemes based on spectral type 
approximations (Canuto et at., 1 987) and is suitable for discretization ol the equations 
considered here. These equations have been given in a number of previous publica- 
tions (Drummond, 1988. Carpenter, 1989; also sec Williams, 1985), and will not be 
repeated here. A special care was taken in generating a fine resolution grid, so that 
the essential features of the reacting compressible field are not concealed by numerical 
errors. 


3 PRESENTATION OF THE RESULTS 

Calculations are performed with different values of the convective Mach number A /, , 
and the heat release parameter Ct\ in order to assess the influence ol these parameters 
on the structure of the layer. In this assessment, all the other non-dimensional 
parameters have been kept constant to isolate the effects of compressibility and 
exothermicity. The simulations were performed on a domain discretized with 
256 x 256 grid points, and some on a 128 x 128 mesh for accuracy checks. A 
constant grid size distribution was employed in the streamwise direction, i.e. constant 
AX, but in the cross stream, the grids were unequally spaced with a heavy mesh 
concentration near the center of the layer (Y — 0) covering the region dominated by 
large scale structures. A schematic of the grid orientation for 128 x 128 simulations 
is shown in Figure 1(b). The routine utilized for grid generation is described in detail 
by Drummond (1988), and will not be discussed here. With the resolution attained, 
it was possible to perform simulations with a Damkohler number of Do, = 100, and 
a Reynolds number of 250 based on the initial vorticity thickness, initial density and 
the free stream velocity. The magnitudes of the Prandt! and the Schmidt numbers 
were assumed to be equal to unity, resulting in an identical value for the Reynolds and 
the Peclet numbers. The results of these simulations, together with a discussion of 
their implications are presented in this section. 

For the purpose of flow' visualization, the plots of the vorticity contours at different 
computational times ( t * = tU x /L x ) are presented in Figure 2. In these plots (and 
those to be presented later), the contours are drawn on a stretched mesh with constant 
A Y, so that the structures and the scales of transport are displayed clearly. This figure 
shows the time development of the vorticity for a case with no heat release at a 
convective Mach number of M i = 0.2. The growth of instability modes in this layer 
is represented by the formation of large scale two-dimensional structures. These 
structures are formed initially in the form of single vortex rollups, followed by the 
pairing of the neighboring vortices to form a large vortex at later times. The results 
portrayed by this figure are similar to those of previous simulations of temporally 
evolving mixing layers (e.g., Riley and Metcalfe, 1980; Givi et ol., 1986). The main 
difference between the results presented here and those of previous simulations is the 
procedure by which the perturbations are added to the mean flow'. The perturbations 
introduced in the present simulations are generated by numerical truncation and 
round-off errors, whereas in previous works explicitly-added harmonic forcing were 
the primary source of perturbations. 

The corresponding plots of the contours of the normalized pressure and the product 
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FIGURE 6 Plots of product mass fraction contours at dilTcrenl time realizations, \f, = 0.8, Cc = 0. (a) 
I* = 6. (b) I* = 8, (c)/* = 10. 


mass fraction for this case are shown, respectively, in Figures 3 and 4. Figure 3 
portrays the pressure response to the formation of large scale structures. As these 
structures are formed, provincial regions of pressure maxima and minima occur. A 
comparison between Figures 2 and 3 indicates that the location of pressure minima 
is at the core of the vortices where the magnitude of the absolute vorticity is the 
highest, and the colony of pressure maxima coincides with the saddle points where the 
absolute vorticity is at a minimum. The formation of large scale structures results in 
a direct enhancement of the combustion product formation, as displayed in Figure 4. 
It is shown that with the development of the large scale vortices, the fluids from the 
two streams are brought into contact, and the chemical reaction occurs at the mixing 
zone with a maximum magnitude of the product mass fraction at the core of the 
vortices. A comparison between Figures 2 and 4 reveals the similarity ol the transport 
of the vorticity and the product mass fraction. This is to be expected since in this case, 
enhanced mixing (caused by the convolution of coherent vortices) corresponds directly 
to an increased combustion and thus an increased product formation. 

The roles played by the compressibility on reaction conversion rate can be quantita- 
tively assessed by examining the magnitude of the product mass traction. The plots 
of the contours of this fraction for different values of the convective Mach numbers 
(keeping heat release rate at Ce = 0) are shown on Figures 5-7 at different time 
realizations. A glimpse at these figures suggests a reverse relation between the mag- 
nitude of the convective Mach number and the extent of large scale mixing and 
subsequent product formation. As A/ t increases, it takes longer for the background 
perturbations to grow' and the layer becomes more sluggish in responding to such 
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FIGURE 7 Plots of product mass fraction contours at different time realizations. A/ = 1 .2. CV - 
i* = 6, (b) r* = 8. (c) t* = 10. 


perturbations. The trend is enhanced as the Mach number is increased and at the 
largest Mach number employed, the rate of growth of the layer and the amount of 
products formed are the smallest. 

The response of the shear layer to increased compressibility can be further ap- 
praised by examining the statistical and the integral properties of some appropriate 
parameters which characterize the flow. For this purpose, the profiles of some of the 
key parameters which include global information about the flow are presented in 
Figures 8-9. In Figure 8, the profiles of the normalized vorticity thickness are 
presented versus the normalized time for different values of the convective Mach 
number. The general trend observed in this figure indicates the reduction of the rate 
of growth of the layer as the Mach number is increased. This is to be expected, in v iew 
of the contour plots shown earlier; increased compressibility infers reduced mixing 
and reduced broadening of generated vorticity. The same is also true in the time 
variation of the normalized total product mass fraction, as evidenced in Figure 9. 

A statistical examination of the fluid mechanical variables at selected computation- 
al times is presented in Figures 10-11. These figures present the cross stream varia- 
tions of the mean and the mean square of the streamwise velocity. The knowledge of 
these quantities is relevant for turbulence closure, and it is the prediction of these 
quantities that is of main interest to turbulence modelers*. With the approximation 
of temporal evolution, the flow is homogeneous in the spatial direction .V, and the 


♦Even though these simulations deal with a purely two-dimensional flow without any three-dimensional 
transport essential for turbulence analysis. 
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FIGURE 8 Normalized vorticity thickness versus normalized time for different values of the convective 
Mach number. 


statistical information is obtained by sampling the data at the grid points in this 
direction. In this way, the ensemble mean and the mean square of a transport variable 
(say <p) are obtained by: 


<<my» = (6) 

<0' : (F)> - ^t(0,-<0» 2 . (7) 

Here < ) denotes ensemble average, and the subscript / and the parameter N 
indicate, respectively, the grid index and the total number of grid points in the X 
direction. 

The profiles of normalized streamwise component of the average velocity are shown 
in Figure 10 at different times. The most significant feature displayed in this figure 
is the steepness of the mean velocity profiles at high Mach numbers. Again, in view 
of the contour plots of the vorticity and the product mass fraction this is to be 
expected, and the increase in the velocity steepness (caused by the reduced growth 
rate), implies lesser rate of mixing and, thus, decreased product formations. This trend 
can also be described by examining the profiles of the mean square of the same 
variables as shown in Figure 1 1 for t* = 6. Note the double hump characteristics of 
the mean square velocity profile at low Mach numbers, with the local maxima (at the 
humps) corresponding to the regions of maximum velocity gradients. Also, note that 
as the magnitude of the convective Mach number is increased, the amplitude of the 
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FIGURE 9 Normalized total product mass fraction versus normalized time for different values of the 
convective Mach number. 


fluctuations decreases, and at M t = 0.8 and M = 1.2, this amplitude becomes very 
small. 

Another interesting characteristic of the increased compressibility is captured by 
examining the plots of pressure contours at high convective Mach numbers. A time 
sequence of instantaneous pressure, portraying the evolution of the layer for convec- 
tive Mach numbers of 0.4, 0.8, and 1.2 are depicted in Figures 12-14. The pressure 
response in Fig. 12 is very similar to that presented for M c — 0.2 case (Figure 3), 
indicating the regions of pressure maxima and minima at the braids and the cores of 
the vortices. At higher convective Mach numbers, however, it is observed that the 
increased compressibility results in steepness of the gradients of instantaneous pres- 
sure and the formation of “eddy shocklets”. These shocklets are initiated at the shear 
zone of the layer, and extend to the outer region of the flow near the boundaries. A 
ratiocination for the formation of these shocklets is provided by noting the increased 
compressibility within the domain in high convective Mach numbers. In these cases, 
the layer is dominated by regions of supersonic and subsonic flows, and in order for 
the flow to adapt to high pressures at the braids, it must go through a shocklet to make 
the proper adjustment. Furthermore, it is observed in these figures that as the 
structures are formed, the regions of high compressibility tend to push the shocklets 
backward. This is evidenced by the pressure contours at small time intervals in Figures 
13-14, which show that the shocklets tend to rotate in an opposite direction to the 
large scale vortical structures. Finally, note that the currents do not necessarily have 
to be supersonic at the free streams, and compression occurs within the flow as a result 
of the formation of large scale structures. This point is more clearly demonstrated by 
examining the contour plots of the instantaneous Mach numbers in Figures 15-16 
(corresponding to some of the cases presented in Figures 13-14). It can be observed 
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FIGURE 10 Profiles of normalized mean velocity (Uy/U t versus F/<5 J 0 for different values of the 
convective Mach number, (a) t* = 6, (b) l* = 8. 

from these figures that for the case of M t — 0.8 (i.e. subsonic free streams) the flow 
at the interior is dominated with localized regions of supersonic ( Ma > 1), and 
subsonic ( Mu < 1) flows. These adjustment from supersonic to subsonic conditions 
is provided by the formation of eddy shocklets. The strength of these shocklets 
become stronger as the convective Mach number is increased (i.e. as the effects of 
compressibility become more pronounced). 

The results of our simulations are consistent with those of experimental measure- 
ments (e.g.. Samimy and Elliot. 1990) in that as the compressibility increases, the 


ORIGINAL PAGE IS 
OF POOR QUALITY 


f IK .11 SIH I n RI \(TI\(| \1I\IN( i 1 \YI R 


4 ') 


(a) 



-50.0 -30.0 -10.0 Y 10.0 300 500 

<^u/|n 


FIGURE 1 1 Profiles of normalized mean square velocity versus YI &,., | u for different values of 

the convective Mach number, (a) t* = 6, (b) i* = 8. 

magnitudes of the turbulence fluctuations decrease. However, the eddy shocklets have 
not been observed in any of the experiments on compressible shear layers, to date. 
This could be a consequence of computational limitations in that the simulations are 
restrictively two-dimensional, whereas the laboratory flow is dominated by small scale 
three-dimensional transport. Whether or not these shocklets are formed in a flow 
under the influence of such transport is to be determined in future three-dimensional 
simulations. 

The effects of heat release can be also assessed by following the same procedure as 
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outlined above. In the following discussions, results arc presenled ol simulations for 
constant rale kinetics with ,\/< = 0.2 and heal release values of Cc = 0. C'c = 1.5, 
and (V = 6; and one set of Arrhenius simulations with Ce = 1.5 and Zc = 10. The 
results of these simulations are portrayed by the plots ol the vonicity thickness versus 
the normalized time for all four cases (Figure 17). This figure shows that the onset ol 
formation of large scale structures is quickest for the non-heat releasing mixing la^er 
and as the effects of chemistry become more pronounced, the response of the layer 
becomes more sluggish. For example note that for the Ce = 0 case, the la^er responds 
to instabilities fairly rapidly, and after the initial thickening there is a jump in the 
magnitude of the thickness at I* ^ 3. This jump indicates the formation of vortical 
structure as shown previously in Figure 2. An increase of the magnitude of the heat 
release (Cc = 1 .5). results in a delay of the vortex rollup, and the jump in the vorlicily 
thickness does not occur until t* ^ 7. Further increase in the magnitude of the heat 
release results in additional delays, and as can be seen for the cases of Cc = 6 and 
Ce = 1.5 (with Arrhenius model), the effects of instabilities do not become apparent 
at all. In these cases, the effect of exothermicity is most pronounced; vorticity rollup 
does not occur and the only cause of vorticity thickening is due to molecular diffusion. 
The contour plots of the vorticity in both of these cases form parallel lines (not show n 
here) which indicate the lack of formation of vortical structures. At larger times, the 
layer becomes too "thick” to respond to perturbations within the flow, and proceed- 
ing in time does not produce any substantial enhancement in mixing except, of course, 
that facilitated by diffusion. The main reason for the severe retardation of mixing for 
the Arrhenius case in comparison with that of constant kinetics model is due to the 
kinetics model employed in the simulations. With the application of the Arrhenius 
reaction model, the rate of increase of temperature is substantially more than that of 
constant rate kinetics, even though the magnitude of the heat release parameter (as 
defined here) is kept fixed. This increase in the local (and global) magnitude ol the 
temperature postpones the growth of the instabilities, and for the case of Zc = 10 (as 
employed here), the layer does not accommodate for the growth of the instability 
modes to form large scale coherent vortices. 

The influence of the heat release on the compositional structure of the flame is 
displayed by examining the amount of normalized total product mass fraction of the 
layer as a function of time. This is shown in Figure 18. An examination of this figure 
shows the influence of heat release on all stages of the development of the layer. At 
initial times, the effect of heat release is a somewhat enhanced product formation, 
whereas at intermediate and final stages a reverse scenario holds. At the initial stages, 
the effect of heat release is to expand the fluid at the cores of the layer. Therefore, it 
is expected to have a larger mixing zone and thus a slightly larger amount ol products 
formed. Flowever, as the heat release increases and the layer thickens, the rate of 
growth of instability modes becomes subdued, postponing the rate ot formation of 
large scale vortices. After the initial stages, the non-heat releasing simulations predict 
a sharp increase in the product formation and as the magnitude of the heat release is 
increased, the time at which such structures are formed gets delayed. The lowest rate 
of product formation is for Ce = 6 simulations in which, as discussed before, the only 
mechanism of mixing is through diffusion. This mechanism of reduced product 
formation is also evidenced by a comparison between the contour plots of the product 
mass fraction with heat release (Figure 19), and those without heat release (Figure 4). 

Further influences of heat release become evident by examining its effect on 
statistical quantities. In Figure 20 the normalized profiles of mean streamwise velocity 
component are presented at different times. This figure shows that at early times, the 
local volume expansion of the fluid caused by exothermicity yields a less steep gradient 
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for the higher heat release rates, whereas at final times increased heat liberation results 
in a relatively more steep gradient of the velocity and, therefore, less mixing. This has 
a substantial influence on the two-dimensional turbulence transport, as indicated by 
the cross stream variations of the mean square velocity presented in Figure 21. As 
shown in this figure, as the flow develops the magnitude of fluctuation is maximum 
for the non-heat releasing mixing layer and, as the exothermicity becomes dominant, 
the amplitude of the fluctuation decreases. For the most significant heat release cases 
(Ce = 6 and the Arrhenius model), the amplitude of the mean square velocity is very 
close to zero indicating virtually no turbulence fluctuations. 

In view of some of the earlier works on DNS of low speed temporally evolving 
reacting mixing layers, there are some points that need to be discussed. These points 
are related, in particular, to the work of McMurtry et al. (1989) and Givi et al. (1986). 
In the former, the simulations of a variable density flow with constant kinetics were 
undertaken, w hereas in the latter, an incompressible layer under the influence of an 
Arrhenius reaction was the subject of main consideration. In the present contribution, 
the response of the layer to heat release for the case of constant rate kinetics is 
consistent with the findings of McMurtry et al. (1989). The major difference between 
present calculations and those of McMurtry et al. is the procedure by which the 
hydrodynamic equations are treated. In this previous work, a low Mach number 
approximation was made, and all the transport equations were filtered to remove the 
acoustic field. This approximation does not allow for simulating flows with high Mach 
numbers and/or high heat release rates. Nevertheless, the general trend of the results 
suggest mixing suppression at elevated heat release rates. It is encouraging to note that 
similar trends are also observed in full compressible simulations in the present 
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contribution. However, in comparison with calculations of Ciivi cl at. (1986). the 
present simulations do not allow the examination of the effects of high stretch rates. 
This is primarily due to accounting for the effects of variable density in present 
simulations which were neglected in the previous work. As displayed by the results, 
the heat release and the associated density variations result in a substantial modification 
of the mixing procedure. In the range of parameters considered, the layer does not 
form large scale structures and, therefore, the magnitude of the stretch rate does not 
become high enough to quench the flame. This indicates the importance of the density 
field in the analysis of flame quenching (Buckmaster cl <//., 1986) which was ignored 
previously (also see the useful comments of Boris in Givi cl at. ( 1986) in this regard). 
In order to investigate the extinction phenomenon, future Arrhenius simulations are 
required in which the effects of heat release are somewhat less dominant. This can be 
facilitated by simulating flows with lower values of the heat release and for various 
magnitudes of the Damkohler and the Zeldovich numbers. Some various combina- 
tions of these parameters are required in future simulations to address this important 
aspect of the non-equilibrium chemistry in detail. This issue is the subject of current 
investigation. 


4 SUMMARY AND CONCLUSIONS 

A high order finite difference algorithm is employed for direct numerical simulations 
of a high speed reacting temporally growing unforced mixing layer. Within the 
framework of temporal approximation, periodic boundary conditions are employed 
in the streamwise direction eliminating the need for providing boundary condition in 
this direction of the flow. With such an approximation, the effects of large scale 
structures and the subsequent coupling between mixing and reaction in the unsteady 
compressible flow are simulated without addressing the asymmetric mixing phenome- 
na which are observed in typical laboratory shear layers. With the accuracy of the 
numerical algorithm (second order in time, and fourth order in space), together with 
fine resolutions (both temporally and spatially), it was possible to perform accurate 
simulations of the reacting field with moderate values of the Reynolds, Peclet, 
Damkohler, Zeldovich, heat of reaction, and the convective Mach number. A full 
compressible code was implemented in order to allow for a reasonable and accurate 
simulation of the high speed reacting field. The main objective of the simulations was 
to examine the isolated effects of compressibility and the heat generated by the 
chemical reactions. A simplified reaction mechanism with idealized kinetics schemes 
was employed. Namely a finite-rate reaction of the type A + B -» Products + Heat 
was the main subject under investigation. Two types of kinetics mechanisms were 
employed, one with constant rate kinetics, and another with an Arrhenius model. In 
both cases, the magnitudes of the Reynolds, Prandtl, Schmidt and Damkohler num- 
bers were kept fixed, but the values of the convective Mach number and the heat 
release parameter were varied. In the constant rate reaction case, the Damkohler 
number is the only parameter which describes the chemistry, whereas in the Arrhenius 
model the Zeldovich number also is an important non-dimensional quantity to be 
considered. 

The results of numerical experiments allowed a reasonable assessment of the roles 
played by the compressibility and the heat release. It was shown that both these 
factors have significant influences on the outcome of mixing, and modify the pro- 
cedures by which the reacting species mix and convert to product species. This, in 
turn, influences the hydrodynamic transport within the flow and modifies the process 
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by which the unsteady large scale vortical structures develop. These influences are 
captured in the numerical simulations by a careful examination of the global struc- 
tures, as well as various integral and statistical properties. 

In the range of parameters considered, it is concluded that heat release and 
compressibility generally result in decreased mixing with the exception of a slight 
increase in mixing caused by volumetric expansion at high heat release rates. This is 
noted by the response of the layer to increased heat release and compressibility. As 
the magnitudes of the heat release parameter and the convective Mach numbers are 
increased, the layer becomes less responsive to background perturbations and it takes 
longer for the instability perturbations to grow. This is evidenced by both flow 
visualization, and by examining the statistical characteristics of the flow. The contour 
plots of various fluid mechanical and chemical quantities indicate that at low heat 
release rates, the instability modes grow at a faster rate, and the large scale coherent 
vortical structures are formed quicker. As the magnitude of the convective Mach 
number is increased, the layer becomes less sensitive to such perturbations. As a 
result, the layer grows at a slower rate and less mixing occurs. This is quantitatively 
demonstrated by examining the temporal variation of the vorticity thickness at 
different convective numbers which present a reasonable assessment of the effects of 
compressibility on mixing. Another feature of the increased compressibility is the 
formation of eddy shocklets as a consequence of increased compression after the 
formadon of large scale structures (even though these structures are not very distinct 
at high convective Mach numbers). These shocklets are formed for both subsonic and 
supersonic free stream flows. As the layer is compressed, the flow in the low pressure 
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regions outside of the mixing core goes through a shock in order to adapt to the high 
pressure regions at saddle points at the braids of the vortical structures. These 
shocklets rotate in an opposite direction to the motion of large scale structures. 

The heat release influences the formation of large scale structures and the subse- 
quent evolution of the layer in two different ways. At the initial stages of mixing, the 
effect of heat release is to increase the volumetric expansion of the fluid. This results 
in thickening of the layer and in increasing the amount of products formed. After the 
initial stage, however, the heat release has a reverse influence. As the magnitude of 
heat release is increased, the layer becomes less responsive to the growth of perturba- 


ORIGINAL PAGE IS 
OF POOR QUALITY 


MK ill SPI tT> RF\< TIN(i \||\l\(. I \YI R 


(> ’ 


1.0 


: U > 

t? 


o.o 


- 1.0 


- 2.0 



-200 


- 10.0 


0.0 


V 


10.0 


20.0 



FIGURE 20 Profiles of normalized mean velocity < U}IU , versus >7 <>, ,,| 0 for different values of the heat 
release parameter, (a) t* = 3, (b) /* = 6, (c) t* = 8. (d) /* = 10. 

tion, and it takes longer for the instability modes to grow. This delays the formation 
of coherent vortices at high heat release and yields a substantial decrease in the 
magnitude of the total products formed. The contour plots of fluid mechanical/chenii- 
cal quantities provide a valuable flow visualization tool in displaying the effect of heat 
release on the formation (or lack of formation) of large scale structures. The cross 
stream variations of the mean and the rms of the same quantities also show the 
augmentation or depression of mixing as a function of the heat release. 

Another important feature of the results of reacting flow simulations is the effect 
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FIGURE 21 Profiles of normalized mean square velocity < U' 2 >/t '] versus >7<VJo f°r different values of 
the heat release parameter, (a) (* = 6. (b) /* = 8, (c) t* = 10. 
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of kinetics modeling. These results indicate t lint lor t lie same values ol the Damkohler 
and heat release parameter, the Arrhenius model ol the reaction conversion has a 
more pronounced effect on mixing hindrance. The exponential temperature depen- 
dent reaction rale increases the magnitude of the temperature much faster, and 
therefore enhances the effects ol heat liberation. In the case with a Zeldovich number 
of 10. is shown that the temperature increase at the initial stages of the growth 
results in a rapid thickening of the layer. At the subsequent stages, the layer does not 
allow for the growth of instabilities, and large scale coherent structures are never 
formed within the time range of the simulations. This results in a substantial decrease 
in the amount of total product formation, and therefore yields a much smaller 
combustion efficiency. Future simulations with different magnitudes of the Zeldoxich 
numbers, as well as the Damkohler numbers and the enthalpy of formation are 
required to examine the effects of these parameters on non-equilibrium behavior ol 
the layer. 

Much remains to be done in future efforts on DNS of high speed reacting flows. 
Several recent reviews are available suggesting the future paths and directions to be 
taken in this important area of research (Drummond, 1991; Givi and Riley, 1991). 
Some of our ongoing efforts in the analysis of reacting high speed shear layers include 
the simulations of three-dimensional flows, investigation ol non-equilibrium effects 
and implementations of more complex kinetics mechanisms. In these future efforts, 
it is planned to examine the effects of secondary instabilities with and without the 
influences of heat release for different convective Mach number regimes. It is also 
intended to examine the response of the layer at different values of the Damkohler and 
Zeldovich numbers in order to make a complementary contribution to the investiga- 
tions done previously for low speed reacting flows (Givi et ai, 1986, Givi and Jou. 
1988). Finally, it is planned to include the effects of more complex kinetics mechan- 
isms on the structure of the flow. In these continuation efforts, the influences ol 
multiple Damkohler and Zeldovich numbers on the formation of large scale struc- 
tures and on small scale three-dimensional random motions are expected to be 
captured by DNS. Some of our other ongoing efforts are concentrated on the 
development and implementation of more accurate numerical algorithms for both 
temporal and spatial discretizations of the transport equations. Two classes of num- 
erical approximations are presently under careful consideration. One is based on 
compact differencing algorithms of Pade type approximations and the other is based 
on the spectral-element methods. Both these methods have proven useful for simulat- 
ing non-reacting and reacting (Carpenter, 1990) and low speed (Givi and Jou, 1988) 
flows. If such methods prove to be expressively more accurate for simulating reacting 
high speed flows, without causing significant inconveniences, they will also be em- 
ployed for future calculations. 
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Abstract 

The effects of compressibility, chemical reaction exother- 
micity and non-equilibrium chemical modeling in a com- 
busting plane mixing layer were investigated by means 
of two-dimensional model free numerical simulations. It 
was shown that increased compressibility generally had a 
stabilizing effect, resulting in reduced mixing and chem- 
ical reaction conversion rate. The appearance of “eddy 
shocklets” in the flow was observed at high convective 
Mach numbers. Reaction exothermicity was found to en- 
hance mixing at the initial stages of the layer’s growth, 
but had a stabilizing effect at later times. Calculations 
were performed for a constant rate chemical rate kinet- 
ics model and an Arrhenius type kinetics prototype. The 
Arrhenius model was found to cause a greater tempera- 
ture increase due to reaction than the constant kinetics 
model. This had the same stabilizing effect as increas- 
ing the exothermicity of the reaction. Localized flame 
quenching was also observed when the Zeldovich num- 
ber was relatively large. 

Nomenclature 

Aj : Arrhenius law constant 

c v : specific heat at constant volume 

C : species concentration 

C t : heat release parameter 

Da : Damkohler number 

E : activation energy 

AH 0 : heat of reaction 

Kj : reaction rate constant 

/o : physical length of computational domain 

M c ‘ convective Mach number 

R : universal gas constant 

t : time 

t m : normalized time 

T : temperature 

U : streamwise velocity 

ibi : species production rate of species : 

Ze : Zeldovich number 


"Research Assistant, Student Member AIAA. 


8 W : vorticity thickness 
Subscripts: 

A } B : chemical species 

x : in x coordinate direction 

y : in y coordinate direction 

co : free stream 


1 Introduction 

Research on the subject of compressible reacting mixing 
layers has been of high priority in recent years. Much of 
this effort has been devoted to the development of high 
speed air breathing flight vehicles. This type of vehicle 
would, according to current proposals, use supersonic 
combustion ramjet (scramjet) engines for propulsion. In 
such an engine, fuel is injected into a high speed airflow. 
The mechanisms of mixing and combustion of this non- 
premixed, high speed, compressible, reacting flow is of 
great complexity. 1 ' 2 

The purpose of this study was to investigate the 
coupling between the hydrodynamic and chemical phe- 
nomena that are present in compressible, reacting mix- 
ing layers. Although experimental studies have provided 
useful results, this avenue of research is often difficult 
and expensive. Numerical simulation has become widely 
accepted as an alternative and supplement to experimen- 
tal investigation. Model free simulations have proven es- 
pecially useful because there is no need to introduce the 
potential inaccuracies of a turbulence model. An exten- 
sive review of current numerical simulation techniques 
and results for reacting flows is provided by Givi. 3 

1.1 Previous Research 

Brown and Roshko 4 found that the turbulent mixing 
layer is dominated by large scale coherent structures, or 
vortices. These structures convect at a nearly constant 
speed and tend to coalesce with neighboring vortices. 
Papamoschou and Roshko 5 continued these experiments 
to determine the effect of compressibility on the spread- 
ing rate of a supersonic mixing layer. They found that it 


1 



is useful l.o study the flow in reference frame Mini, trav- 
els with the flow at the same speed as an average large 
scale structure. A parameter which quantifies the com- 
pressibility in the flow was proposed as the convective 
Mach number, M c \ defined as the Mach number of the 
flow with respect to the above mentioned frame of ref- 
erence. A direct correlation was found between M c and 
the stability of the flow. Lele 6 replicated the results of 
Bapamoschou and Roshko 5 by means of direct numerical 
simulation of a two-dimensional mixing layer. In addi- 
tion, he noted the appearance of eddy shocklets within 
the flow for M c >0.7. 

McMurtry, et a/. 7 studied the effects of an exother- 
mic chemical reaction on the shear layer via direct nu- 
merical simulation. In this study, an approximate set of 
equations that are valid for low Mach number flows was 
used. It was found that the heat liberated from a chem- 
ical reaction causes the layer to grow at a slower rate 
that of on-heat releasing flow. These results agree with 
those obtained experimentally. In addition, direct nu- 
merical simulations concerned with compressibility and 
heat release effects on a compressible mixing layer were 
performed by Givi et a/. 8 

The phenomenon of flame extinction in non- 
premixed flames has been the topic of theoretical and 
experimental study, although direct numerical simula- 
tions of such flows have been somewhat limited. Recent 
reviews of some of the prevalent theories regarding the 
structure of turbulent non-premixed flames, as well as 
some of the experimental and numerical work in that 
Held are provided by Bilger 9 and Peters. 10 

1.2 Scope of Present Research 

The scope of the present work is to examine the effects 
of compressibility ami chemical reaction exotherm icit.y 
on a reacting plane mixing layer. An examination is 
also made of the non-equilibrium effects of the chemical 
kinetics on the structure of a flame. These are accom- 
plished by direct numerical simulation of an unsteady 
two-dimensional shear layer. The governing equations 
are integrated via high order finite difference methods, 
and no turbulence modeling is employed. Since the focus 
of this study is the general physics of this type of flow, 
physical modeling is kept as simple as possible so that 
the above mentioned hydrodynamic and chemical effects 
may be isolated. The use of simplifying assumptions has 
the added advantage of saving considerable amounts of 
computational resources, without reducing the accuracy 
or validity of the results. 


2 Method of Investigation 

A two-dimensional version of the SPARK computer 
code 11,12 was used to investigate a chemically react- 
ing, compressible, temporally developing mixing layer. 


The full Navier-Stokes equations i\s well as the applica- 
ble chemical species conservation equations (see 11 or ia 
for more detail) were solved via numerical integration. 
The methods used were the Gottiicb-Turkel finite dif- 
ference scheme 11 and a dissipative compact parameter 
method. 15 Both of these algorithms are of fourth order 
spatial accuracy and second order temporal accuracy. 

One of the primary assumptions made in these sim- 
ulations is that the mixing layer is temporally develop- 
ing. That is, the reference frame of the simulations is 
defined to be moving with a velocity equal to that of 
the mean flow. The advantages of this approximation 
are twofold. First, with the temporal assumption the 
inflow and outflow boundary conditions can be assumed 
to be periodic. The eliminates the problem of specifying 
downstream boundary conditions. Second, the temporal 
assumption means that only a relatively small region of 
the flow is being simulated. This region is then followed 
in a Lagrangian sense as time progresses. This results in 
considerable savings in computational resources, which 
means that the flow may be simulated in greater detail. 

The chemical reaction in the flow is assumed to be 
of a simple, irreversible, second order type of the form 

A + B — ► Product + Heat (1) 


The reaction is characterized by the kinetics mechanism, 
which is given by the single step model of 


w^I< s C A C B ( 2 ) 

where Ca and Cq represent the concentrations of the re- 
acting species and are assumed equal at the free streams, 
i.e. Ca^ = C Boo = Coo- Kj is the reaction rate con- 
stant, and can be normalized to form the definition of 
the Damkohier number, Da: 


Da = 


HjCqo 

Hoo / S w \q 


(3) 


In the present study two types of chemistry models were 
used; constant rate kinetics (i.e. constant Kj ) and an 
Arrhenius type model in which Kj varies with the tem- 
perature. This is written as 


Kj — A jt ^T^oo 


(4) 


where Aj is the pre-exponential factor and Ze is the 
Zeldovich number, defined as 


Ze = 


E 

KToo 


(5) 


Here E is the activation energy and R is the universal 
gas constant. Too is the free stream temperature. When 
the Arrhenius kinetic model is used, the pre-exponential 
factor Aj replaces Kj in the definition of Da (Eq. 3). 

Combustion exotherm icity is measured by the en- 
ergy liberated by the chemical reaction, AH 0 , The 
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Figure 1: Schematic diagram of a temporally evolving 
mixing layer. 


magnitude of this energy is parameterized by a non- 
dimensional heat release parameter Ce, defined by: 


-A H° 


Thus, Ce = 0 corresponds to a non-heat releasing chem- 
ical reaction. 

The flow was initialized with a hyperbolic tangent 
streamwise velocity profile, with reactant A on the top 
half of the layer with a free stream velocity of -\-Uoo and 
reactant B on the bottom half with a free stream veloc- 
ity of — Uoo- This flow configuration is shown schemat- 
ically in Fig. 1. There was no initial fluid motion in 
the transverse direction, and the pressure was assumed 
to be initially constant throughout the flowfield. De- 
pending on the problem simulated, the temperature was 
either assumed to be initially constant or had an initial 
Gaussian distribution. In most cases, the mixing layer 
was perturbed by numerical truncation errors. However, 
when necessary, harmonic forcing was explicitly added 
to trigger the flow instabilities. For a majority of the 
simulations, a 128 x 128 point grid was used, with grid 
compression in the areas of the flow with the sharpest 
gradients. In those cases where the gradients were too 
strong for the details of the flow to be resolved, a 256 x 
256 point grid was used. 

All parameters are normalized when appropriate by 
initial or free stream conditions. Time is normalized by 


3.1 Effects of Compressibility 

The efTect of compressibility on the mixing layer was in- 
vestigated by varying M C) the convective Mach number 
over a range of values from M c — 0.2 to M c = 1.2. The 
vorticity thickness vs. normalized time for different val- 
ues of M c is given in Fig. 2. This figure clearly shows that 
the growth rate of the vorticity thickness is decreased 
with increased compressibility. Compressibility also af- 
fects the time needed for the layer to roll up into a vortex, 
which is shown in the figure as a jump in vorticity thick- 
ness. Valuable information may be obtained from exam- 
ination of the variation in the transverse direction of the 
normalized average streamwise component of the veloc- 
ity and the normalized mean square streamwise velocity 
fluctuations. The average streamwise velocity profile is 
shown in Fig. 3. The most significant feature of these 
curves is the steepness of the mean velocity profiles at 
high convective Mach numbers. This shows a sharper 
velocity gradient across the layer, implying a lesser rate 
of mixing. The degradation of mixing has the effect of 
retarding the progress of the chemical reaction. The sup- 
pression of turbulence fluctuations with compressibility 
is shown by the profiles of the mean square of the fluc- 
tuating velocity in Fig. 4. A marked decrease in the 
amplitude of the fluctuations can be observed. 

Another feature of the increased compressibility is 
the formation of eddy shocklets as a result of increased 
compression after the formation of large scale structures. 
These shocklets are formed for both subsonic and super- 
sonic free stream flows in order to adapt to high pressure 
at the braids. An example of this can be seen upon ex- 
amination of Mach number contours, shown in Fig. 5. 
This phenomenon has not been observed in any of the 
experiments on compressible shear layers, to date. This 
might be caused by the fact that the simulations are two- 
dimensional, whereas a laboratory flow is dominated by 
small scale three dimensional transport. It is anticipated 
that future three dimensional simulations will clarify the 
matter. 

3.2 Effects of Reaction Exothermicity 

Calculations were also performed to assess how reaction 
exothermicity aflects the mixing layer. Results are pre- 
sented of simulations for constant rate kinetics with heal 
release values of C e = 0, 1.5, and 6. One set of calcula- 
tions was performed using Arrhenius type kinetics with 
C t = 1.5 and Ze = 10. 

The results of these simulations presented in the 
form of plots of the vorticity thickness versus normalized 
time (Fig. 6). This figure shows that the rate of growth 
is highest for the no-heat release case ( C e = 0), and, as 
the heat release parameter is increased, the growth of 
the layer slows. The relatively smooth regions of vor- 
ticity thickness growtli may be attributed to diffusion 
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Figure 2: Normalized vorticity thickness versus normal- Figure 4: Profiles of normalized mean squared stream- 
ed time for different values of the convective Mach wise velocity for different values of the convective Mach 

number. number at time t* = 8. 


£LO 


Figure 3: Profiles of normalized mean streamwise veloc- 
ity for different values of the convective Mach number at 
time t* = 8. 


Figure 5: Plot of Mach number contours at time t 9 

11 . 66 , M c = 1 . 2 . 
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Figure 6: Normalized vorticity thickness versus normal- 
ized time for different values of the heat release param- 
eter. 

thickening, and a jump in vorticity thickness represents 
vortex roll-up. For the C t — 0 case, the layer responds 
to background perturbations fairly quickly, and vortical 
structures are formed at /* % 3. An increase in the 
magnitude of the heat release results in a delay of vor- 
tex roll-up, and the jump in vorticity thickness does not 
occur until t* % 7. Further increase of the heat release 
parameter results in additional delays, as can be seen for 
the case of C e = 6. This behavior is also observed for 
the Arrhenius model with C e = 1.5. In these two cases, 
the effects of exotherm icily are most pronounced; vor- 
ticity roll-up does not occur at all, and the only growth 
in the thickness of the mixing layer is due to molecular 
diffusion. 

The influence of the heat release on the structure 
of the flame is demonstrated by examining the product 
thickness of the layer, defined by the normalized total 
product concentration of the layer as a function of time 
(Fig. 7). At initial times, the effect of heat release is 
an enhanced product formation, while the reverse ap- 
plies at the intermediate and final stages. Initially, the 
efTect of heat release is to expand the fluid at the cores 
of the layer. A large mixing zone is formed, which re- 
sults in an enhanced reaction and an increased product 
formation. This explains the increased initial product 
thickness. However, as the heat release increases and 
the layer thickens, the growth of the instability modes 
become subdued, postponing the formation of the large 
scale vortices. After initial times, the non-heat release 
(C e = 0) and C e = 1.5 cases predict a sharp increase in 



Figure 7: Normalized product thickness versus normal- 
ized time for various values of the heat release parameter. 

the product thickness. This is caused by the dynamics 
of the large scale structure formation. The mechanism 
if roll-up causes the reactants to be entrained into the 
layer from their prospective free streams. This produces 
an increase in the extent of mixing, reaction, and prod- 
uct formation. The lack of roll-up in the C e = 6 and 
Arrhenius cases means the only mixing process is molec- 
ular diffusion. Since in these simulations the diffusion 
mechanism is not as efficient as the roll-up of vortices in 
the mixing of reactants, the product thickness is corre- 
spondingly lower. 

3.3 Flame Extinction 

Another area of study was the nonequilibrium effects 
leading to local flame extinction. Simulations were run 
for large values of the Zeldovich number for an Arrhenius 
kinetics model and for a constant rate kinetics chemistry 
model as a control. To obtain the vortex roll-up neces- 
sary to this analysis, explicit harmonic forcing was added 
to the initial conditions of the flow. A forcing amplitude 
of 0.5% of the mean velocity was used. The most unsta- 
ble frequencies were found by a linear stability analysis 
of an equivalent incompressible flow. It was found that 
if the flow was excited only by the most unstable fre- 
quency, two vortices form but they do not pair. If the 
first subharmonic is added, the vortices undergo pairing. 

In the case of large Zeldovich numbers, it was ob- 
served from the simulations that, after roll-up and pair- 
ing, the flame undergoes local extinction at the braids. 
This behavior is depicted in contour plots of the reaction 
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Figure 8: Plot of reaction rate contours, Ze = 20, i* = 1. 


rate at times t* = 1, 1.5, 2.25, and 2.5. These contours 
are shown in Fig. 8 for time t* = 1, when the vortices 
have just formed. 

The reaction rate is highest along the mixing sur- 
face of the layer, that is, in the center of the intervortex 
braids. The contours at time t* = 1.5 (Fig. 9) portray 
the behavior at the initial stages of pairing. The reaction 
rate has begun to decrease in the braids of the emerged 
vortex. At time /* = 2.25, shown in Fig. 10, the layer 
has completed the pairing process and strong gradients 
are observed in the vortex braids. Since the braids arc 
"stretched” as the vortex rotates, they are the area of 
the highest strain. The flame at the onset of extinction is 
shown in Fig. 11. Note that the reaction rate has gone to 
zero at the braids. Product concentration contours (not 
shown) indicate that very little product exists in these 
extinguished regions, demonstrating that the flame did 
not quench due to depletion of reactants. At this point 
in time, the flame is not continuous, forming what will 
be called for lack of a better term, a "flame eddy.” 

This extinction phenomenon has been explained by 
Peters 10 as follows: At the regions of high strain, the 
reactants are supplied at a faster rate than they can be 
consumed by the flame. Thus the local temperature in 
that area drops and the flame becomes very rich with 
reactants. As a result, the flame is quenched in that 
area. If a fast chemistry model or an equilibrium chem- 
istry model is used, the extinction mechanism can not 
be captured. Investigation of such phenomena requires 
finite-rate chemistry simulation such as that presented 
here. 



Figure 9: Plot of reaction rate contours, Ze — 20, t* = 


1.5. 



Figure 10: Plot of reaction rate contours, Ze = 20, t m = 
2.25. 
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Figure 11: Plot of reaction rate contours, Zt = 20, t* = 
2.5. 

4 Conclusions 

Direct numerical simulations of a compressible, tempo- 
rally developing, reacting plane mixing layer were per- 
formed. Several simplifying assumptions were made so 
that the effects of the variation of isolated parameters 
could be studied in detail. In particular, the chem- 
ical reaction was assumed to be of the general form, 
A + D — ► Products + Heat , and thermodynamic proper- 
ties were assumed constant and identical for all species. 

Studies of various flow phenomena were performed 
by varying one representative nondimensionalized pa- 
rameter while keeping all other parameters constant. 
The convective Mach number, M Ct was used to de- 
scribe compressibility, the heat release parameter, C e , 
represented the exothermicity of the reaction, and the 
Damkohler number and the Zeldovich number described 
the rate of the reaction. 

The simulations concerning the effect of compress- 
ibility on the mixing layer showed a direct correlation 
between increased compressibility and increased stability 
(along with reduced turbulence). When the convective 
Mach number was increased, the rate of growth of the 
mixing layer, which was measured by the growth of the 
vorticity thickness, was markedly reduced. Degradation 
of the development of the streamwise fluctuating veloc- 
ity and mean velocity profiles was also noted. For high 
compressibility cases, “eddy shocklets” were observed 
within the flow. This has been previously reported in 
two-dimensional simulations, but has not been observed 
in experiments or in three-dimensional studies to date. 


An expansion of this work into three dimensions is sug- 
gested for future work in this area. 

Increased exothermicity was observed to slow the 
growth of large scale structures. At the initial stages of 
development, high heat release increased the amount of 
product formed via volumetric expansion of the core of 
the layer. However, the heat release caused the layer to 
be Less responsive to perturbations, reducing the growth 
of the layer at later stages. The overall effect of increased 
heat release, therefore, was to stabilize the flow and to 
decrease the extent and efficiency of the reaction. 

The selection of chemical kinetics model was shown 
to have significant effect on the development of the flow. 
The introduction of an Arrhenius chemistry model had 
a stabilizing effect, thus degrading the progress of the 
reaction. Non-equilibrium chemistry effects due to the 
chemistry model was also studied. It was found that at 
high Zeldovich numbers, the flame would be quenched 
at regions with large local values of the strain rate. 
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ABSTRACT 

Direct Numerical Simulations (DNS) are performed to investigate the phenomena of 
mixing in a decaying two-dimensional, homogeneous turbulent flow under non-reacting 
and reacting unpremixed conditions. In the reacting case, a second order stoichiometric 
irreversible reaction of the type A + B — *■ Products is considered. An infinitely fast 
chemical reaction rate is assumed without including the effects of non-equilibrium chem- 
istry. The data obtained by DNS is utilized to examine the behavior of the scalar field 
in a stochastic manner. This involves the extraction of Probability Density Functions 
(PDF’s) of the scalar field from the DNS results. The examination of the data in this 
manner provides a useful mechanism for assessing the range of validity of PDF methods 
for turbulence modeling in Reynolds averaged equations, and for subgrid closures in 
Large Eddy Simulations (LES). 

Simulations are performed of flows with various levels of compressibility to examine 
its effects on the structure of the flow field. In the low compressibility regime, the 
flow exhibits features similar to those observed in previous incompressible simulations, 
whereas in the high compressibility case, shocklets are formed. In both cases, the 
simulated results indicate that the PDF of a conserved Shvab-Zeldovich scalar quantity, 
characterizing the mixing process, evolves from an initial double-delta distribution to 
an asymptotic shape that can be approximated by a Gaussian distribution. During this 
evolution, the PDF can be approximately described by a Beta density and also exhibits 
features compatible with those predicted by the mapping closure of Kraichnan. The 
main difference exhibited between the two cases is the time scale of decay of the rms 
of the scalar by which the PDF is characterized. A discussion of the results is provided 
with a focus on their relevance towards turbulence modeling and LES closure. A brief 
review is also provided of recent related contributions in PDF modeling and LES in 
order to provide a framework for anticipated continuation of the present work. 

1. INTRODUCTION 

In review of recent computational work on turbulent reacting flows, Drummond 
(1991), Oran and Boris (1987) and Givi (1989) indicate the need for further develop- 
ment in both the methodology and the implementation of Direct Numerical Simulations 
(DNS) and Large Eddy Simulations (LES). According to these reviews, the results of 
previous and ongoing works towards the development of advanced numerical algorithms 
for simulating reacting turbulent fields have been successful, thus warranting their con- 
tinued utilization for future implementations. With the broad knowledge gained to 
date, it is now widely accepted that foreseeable developments in advanced computa- 
tional facilities will not be sufficient to relax the restriction of DNS to flows having 
small to moderate variations of the characteristic length and time scales (Rogallo and 
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Moin, 1981; Reynolds, 1990; Schumann and Friedrich, 1986, 1987, Givi, 1989). Nev- 
ertheless, DNS can be (and have been extensively) used to enhance our understanding 
of the physics of chemically reacting turbulent flows by providing specific information 
concerning the detailed structure of the flow, and by furnishing a quantitative basis 
for evaluating the performance of turbulence closures. In this sense, DNS have pro- 
vided an effective tool for such studies and have established useful guidelines for future 
investigations. 

LES appear to provide a good alternative to DNS for computing flows having ranges 
of parameters similar to those encountered in practical systems (Reynolds, 1990; Hus- 
saini et ai, 1990, Ferziger, 1981, 1983, Schumann and Friedrich, 1987). This approach 
has a particular advantage over analogous Reynolds-averaging procedures in that only 
the effects of small-scale turbulence have to be modelled. However, despite their success 
for the analysis of both incompressible (Ferziger, 1981) and compressible (Erlebacher 
et al . , 1987) flows, they have seldom been employed for calculations of reacting turbu- 
lence. This is primarily due to appearance of enormous unclosed terms in the transport 
equation governing the evolution of the filtered quantities in LES. In these equations, 
the correlations amongst the scalar variables must be accurately modelled in order to 
provide a realistic estimate of the reactant’s conversion rates. 

A methodology which has proven useful for predictions of turbulent reacting flows is 
based on the Probability Density Function (PDF) or the joint PDF of the scalar variables 
which characterize the field (O’Brien, 1980; Pope, 1979, 1985, 1990a). A particular 
advantage of the method is the fundamental property of the probability density, namely 
that it provides the complete statistical information about the behavior of the variable. 
This, in comparison with moment methods, provides a unique advantage in that once the 
distribution of the PDF is known, the mean and higher order moments of the reaction 
rate (or any other functions of the scalar field) can be directly determined without a need 
for further closure assumptions. The implementation of PDF methods in conjunction 
with the Reynolds- and/or Favre- averaged form of the equations of turbulent combustion 
has proven very successful (For an excellent recent review see Pope (1990a)). During 
the past 15 years, these methods have been used for simulating a variety of reacting 
flow configurations and, in many cases, the predicted results have been indicative of the 
higher accuracy of these methods in comparison with conventional moment methods. 
Based on this indication, it is safe to state that presently, from the statistical point 
of view, PDF methods are the most attractive choice for the description of chemically 
reacting turbulent flows. 

A similar scenario may be surmised to be the case when the PDF methods are used 
in the context of subgrid closures for LES. In this case, the PDF (or the joint PDF) 
of the scalar variables within the ensemble of grids (identifying the large scales of the 
flow) may uniquely determine the subgrid average value of the reaction rate (or other 
functions of the scalar variables). This PDF, in conjunction with an appropriate model 
for the hydrodynamic fluctuations within the subgrid, may facilitate the implementation 
of LES for the analysis of reacting turbulence. 

Our main goal in this work is to initiate an assessment of the use of PDF methods 
for their possible future implementation in LES, as well as to further our investigation of 
the PDF methods in modeling turbulent combustion problems. The approach followed 
is similar to those previously used in the a priori analyses. Namely, the data obtained 
by DNS are utilized to construct the appropriate PDF’s that can be used for either 
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turbulence modeling or subgrid closures. We are currently in a preliminary stage of our 
work, especially in tasks directed at subgrid modeling. 

In order to establish a framework for the discussion of the results, a review is 
provided of PDF methods and subgrid closures in Section 2. This review is rather 
brief and is presented solely for the purpose of providing a reference for the remainder 
of the chapter. However, the reader is referred to other collective works and survey 
articles which discuss these statistical methods in greater detail. Also, the discussion 
on PDF methods focuses on the treatment of scalar quantities alone, without including 
the treatment of scalar-velocity PDF’s. The specification of the problem under current 
scrutiny is provided in Section 3, in which the relevance of some of the characteristic 
variables are expressed, together with a description of the geometrical configuration 
and the initialization procedure. The results of simulations are presented in Section 4, 
in which we place an emphasis on the use of stochastic methods for modeling of the 
mixing phenomenon. This chapter ends in Section 5, where we summarize our findings, 
and also provide some speculations and suggestions for future work, some of which are 
currently under way. 


2. BACKGROUND REVIEW 


2.1. PDF METHODS: 

In the majority of previous work on the computational treatment of turbulent 
combustion, statistical methods have generally been employed. These methods involve 
the representation of physical variables in a stochastic sense, and in conjunction with 
appropriate transport equations, predict the approximate behavior of the flow field. In 
the framework of such methods, the physical quantities axe treated as random variables, 
and the transport equations describing the evolution of these quantities are used in such 
a way as to yield their stochastic variations. At present, two general types of strategies 
are identified which can be utilized for statistical treatment of turbulent combustion 
problems (Libby and Williams, 1980): (1) the moment method, and (2) the Probability 
Density Function (PDF) method. In the first approach, “averages” (or “means”) of 
the relevant physical variables are introduced, and the governing transport equations 
are subsequently averaged. These averages are generally defined in such a way to yield 
the desired mathematical properties and axe usually introduced in the form of time-, 
ensemble-, or Favre-averages. The second approach, the PDF approach, relies directly 
on the PDF of these variables by which the desired moments can be obtained. 

Moment methods usually involve the determination of the statistical means by 
virtue of solving transport equations for the various unclosed moments of the equations. 
These equations provide, in a sense, a description of the variables at one level higher 
than the equations for the mean values (first moments). This results in explicit transport 
equations for the second order moments. These equations, however, are still not in a 
closed form, due to the appearance of the unclosed third order correlations wjiich require 
further modeling. Regardless of the degree of truncation of the higher moments, the 
closure problem always remains, and transport equations at any level require modeling 
at higher orders. 

The primary advantage of moment methods is their relative ease of use in en- 
gineering applications. Most models suggested in the literature can be implemented 
directly into computer codes developed for laminar flows. This allows the engineer to 
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economically obtain approximations to combustion problems for many different sets of 
flow conditions. The results obtained in this way have, in many cases, exhibited good 
agreement with experimental measurements. This has been a major factor in justifying 
the use of moment methods in turbulent combustion and it appears that these methods 
will remain popular for engineering applications for the near future. 

The closure based on the second approach, the PDF method, has proven very 
useful in the theoretical description of turbulent reactive flows (Hawthorne et al., 1949). 
The philosophy of this approach is to consider the transport of the PDF’s rather than 
the finite moments of the statistical variables. The key advantage of the method is 
the fundamental property of the PDF, namely that the mean reaction rate (and any 
other functions of the scalar) can be evaluated directly from the PDF. Representing the 
scalar field involving a species by $ = (<p l , (f > 2 , . . . , <f> a ), and the reaction conversion rate 
of species a by w Q ($), the mean of this rate (denoted by < >) can be written as: 

/ +oo 

'P 1 (E)w a (E)dH (1) 

-OO 

In this equation, 5 represents the scalar space and 'Pi(E) is the PDF of the scalar 
variable at a point, defined so that: 

Vi(E)dE = Probability (E < $ < E + d=) (2) 

where dZ = (d^ 1 , d£ 2 , . . . , d^) denotes an elemental hypervolume at E. 

The scalar PDF, Pi(E), defined in Eqs. (l)-(2), is in its simplest form in that it 
contains information only about the scalar variables, 3>. Also, it governs the probability 
distribution only at a single point. However, it contains much more information than 
is required to determine the mean value of any functions of the random variables $. 
Therefore, its evaluation in comparison to direct evaluation of < u> Q > (if it could have 
been done) is understandably more difficult. Nevertheless, since Vi{Z) includes all the 
statistical information about the scalars, its determination is in many ways more useful 
than that of the mean values. 

A systematic way of evaluating a PDF is to obtain and solve a transport equation 
governing its evolution. A transport equation for / P\( < z ri i = £ ) can be 

constructed by relating the PDF to $(x,t) in terms of Dirac-delta functions (Pope, 
1979, 1985, 1990a; O’Brien, 1980, 1981, 1986): 

Vi(E) =< g(x,t) > (3) 

where Q is the single-point fine grained density defined by: 

e = n 6 ^ a ~ *“) (4) 

G — 1 

From the appropriate conservation equations for the scalar field and Eqs. (l)-(4) 
one can obtain transport equation describing the evolution of the PDF. Limiting the 
formulation to that of a constant density Fickian diffusion flow, we have (see O’Brien, 
1980 for a complete derivation): 
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f+<V'V«>=-^M) + r<^/|> (5) 

Here, V is the velocity field and T is the diffusion coefficient of the species. From 
this equation, it can be seen that the effects of reactivity and diffusivity (the terms 
on the RHS) are to transport the PDF into the composition space (5), whereas the 
convective transport is a physical space phenomenon. Furthermore, it is shown that the 
influences of w 0 ($) in the composition space appear in closed form. However, models 
Eire needed for the closures of the molecular mixing term (second term on the RHS) and 
the turbulent convection term (second term on LHS). 

Since the term that includes the interaction between the scalar variables through 
the reaction rate u Q appears in a closed form in Eq. (5), in what follows we limit 
the discussion to the transport of a single reacting scalar (i.e. a — 1) which is being 
convected by the velocity field. In this frame, the diffusive transport may be expressed 
as (O’Brien, 1980, 1981): 

< >= vi [£«2|6)Pi(fi)] (6) 

Here, .E^fol^i) represents the expected value of the scalar at x 2 from its known value at 
Xj . This equation explicitly reveals a fundamental difficulty associated with the PDF for- 
mulations, namely that the transport equation for a one-point PDF (V i(£, x, t )) requires 
information about the second point, x 2 . In general, n-point PDF’s, V n (£i, £ 2 , • • • ,£n, 
x l7 x 2 , . . . , x n , t) require information about the n + 1th point, x n+1 . This is further 
demonstrated by considering the transport of V n in the compositional space (Ievlev, 
1970, 1973; Kuo and O’Brien, 1981): 
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C (q) = riz'mx B+1 -x U) V| n+1 [£(£«+ iKi,6,---^n,x 1 ,x 2 ,...,x n )] (8) 

and E(U+ i|£i,£ 2 ,-.-,£n,Xi,x 2 ,... ,x n ) is the expected value of £ n+l at the point x n+1 
based on all the values of £i,£ 2 ,- • • ,Cn at x^Xj, . . . ,x n . The dependence of C (?) on the 
n + 1th point in the transport equation for V n indicates the need for a closure model 
which relates V n +\ to V n . 

Despite the attractiveness of the PDF methods in obtaining a closed representation 
of the reaction term, there is another difficulty with the implementation of the method. 
This difficulty is caused by the increase of the dimensionality of the PDF as the number 
of scalars and/or the number of statistical points increase. This imposes a severe limit 
on the maximum number n in the equation for V n • The majority of previous work on 
PDF’s have closed the equations at the first level (n = 1), and only recently has closure 
at the n = 2 level been attempted. Both these closures are reviewed in detail by Givi 
(1989). 
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Q.2. LES AND SUDGRID CLOSURES: 

Despite the present capability of modern supercomputers in allowing calculations 
with more than one million grid points, the range of length and time scales that can be 
resolved by DNS is substantially smaller than those of turbulent flows of practical inter- 
est. In DNS the largest compj^table scales are limited by the size of the computational 
domain and the “turn-over” time of the large scale structures; the smallest resolvable 
features are limited by the molecular length and the time scales of the viscous dissipation 
and chemical reactions. DNS therefore, in comparison to turbulence modeling, find its 
greatest application in basic research problems in which the scales of the excited modes 
remain within this band of computationally resolved grid sizes. This implies that for an 
accurate simulation, the magnitudes of the viscosity and the diffusivity must be large 
enough to damp out the unresolved scales, and the size of the computational time step 
must be kept small enough to capture the correct temporal evolution of the turbulent 
flame. In this context, since the instantaneous behavior of the flow variables is directly 
simulated at all spatial points at every instant, the computed data can actually be used 
to provide a quantitative basis for evaluating the validity of turbulence models and for 
assessing the performance of subgrid closures in LES. 

The limitations associated with DNS may be alleviated, to an extent, by pre- 
filtering the transport equations (Aldama, 1990), a practice implemented in LES. This 
is effectively equivalent to eliminating the scales smaller than those resolvable within a 
given mesh. In this way, the variables can be represented on the number of grid points 
that are available for the simulations. This filtering procedure facilitates the simulations 
of flows with larger parameter ranges on a coarser grid. The disadvantage is that some 
modeling is required for the closure of the scales excluded by filtering. The selection 
between DNS and LES (and turbulence modeling for that matter) is dependent on the 
type of flow being considered, on the range of physical parameters that characterize the 
turbulent field, and the degree of desired accuracy. 

A straightforward implementation of LES involves the decomposition of the trans- 
port variables into “large scale” and “subgrid scale” components. The former is related 
to the large scale eddies in the turbulent field, whereas the latter is the component con- 
taining the small scale fluctuations. The pre-filtering of the variable $(x, t ) is performed 
by means of the convolution integral (Ferziger, 1977, Aldama, 1990): 

$(x,/) = J J j F£(x -x')$(x',0<&' (9) 

where Ft is an appropriate filter function with a characteristic length A, and the inte- 
gration is over the entire flowfield. The remaining portion of $ from the filtered quantity 
is defined as the subgrid-scale field and is represented by: 


$"(x,/) = $(x,<) - $(x,t) (10) 

The magnitude of pre-filtered averaged $ and its subgrid component are obviously 
dependent on the filter type (function Ft) and its size A (for reviews see Ferziger, 1981, 
1982, 1983). The simplest choice would be a box filter with: 


-*'> = { 0 , 


if |x — x'l < A; 
Otherwise. 


( 11 ) 
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where A = (A r , A y , A*) represents the dimensions of the box, and the length in each 
direction should be larger than the original grid spacing in that direction. Other types 
of filters have also been suggested in the literature, each of which with computational 
and physical advantages as well as limitations (see Aldama, 1990 and Schumann and 
Friedrich, 1986, 1987). 

Implementation of LES as a practical tool involves a combination of DNS for the 
filtered portion of the transport variable $, and modeling of the subgrid-scale portion, 
0". The transport equations for the large scale components are obtained by filtering 
the instantaneous transport equations by means of implementing the filter defined by 
Eq. (9). In the resulting filtered equations, analogous to those in Reynolds averaged, 
the equations representing the large scale field contain unclosed terms involving the 
fluctuation of the subgrid components. The methodology practiced to date in dealing 
with these fluctuations is very similar to that employed in Reynolds averaging proce- 
dures. One may be able to develop and solve transport equations for the moments of 
fluctuations (Ferziger, 1981). These equations, analogously, contain some higher order 
moments which are needed to be modelled. The number of these terms are usually more 
than the corresponding ones that appear in Reynolds averaged transport, because some 
items that are zero in the Reynolds averaging approach are non-zero when filtering is 
used. Correspondingly, the tasks required for modeling the subgrid fluctuations would 
be somewhat more complicated than those in turbulence modeling. Nevertheless, since 
the small scales of turbulence are believed to exhibit a more “universal” character is the 
main reason to anticipate that the approach based on subgrid modeling would be more 
promising than the procedures based on ensemble averaging closures. 

Within the past two decades, there have been many attempts to fine-tune the 
models by optimizing the closure for subgrid fluctuations and the type of filtering to 
identify the large scale components of the variables. In most of these efforts, the closure 
of hydrodynamic fluctuations has been the subject of major concern. There has also 
been some results for the passive scalar simulation and modeling of velocity-scalar cor- 
relations. However, no attempts have been made for the treatment of reacting scalars 
and the treatment of scalar-scalar interactions. This is probably due to the consensus 
among the researchers in this field that until an accurate subgrid model is constructed 
to represent the evolution of non-reacting scalar variables, the extension to reactive flow 
simulations will be a difficult task. 

With the new developments and progress in PDF methods and the advantages of- 
fered by such schemes over moment methods, one may anticipate that these methods 
may also prove useful in LES. In this case, the influences of the chemical reactions 
and the subgrid scalar-scalar correlations can be included by means of considering the 
PDF’s of the fluctuating scalars within the computational grid. The advantage of using 
these methods for the subgrid closure is apparent since once these PDF’s are known, 
any statistical quantity related to scalar fluctuations can be subsequently determined. 

■ This determination, although an ambitious task, is not impossible. The obvious and 
simplest choice is to follow an approach based on guessed PDF methods.* Similar to 
Reynolds averaging, the first two subgrid moments of the scalar variables can be solved 
by LES and then the shape of the PDF can be parameterized based on these two mo- 
ments. This parameterized approximate distribution can be appraised by a comparison 
with the “exact” shape constructed via DNS results. Similar to common practice in 
a priori assessments (Erlebacher et a/., 1987), the PDF distribution can be specified 



8 


by performing two sets of calculations; one with a coarse mesh, the other with a fine 
mesh. The results obtained from the fine mesh simulations can be used to construct 
the PDF distribution which in turn may be used in extensive subsequent simulations 
on the coarse mesh. In this setting, calculations with large transport parameters which 
otherwise could not be simulated on the coarse grids, are possible. 

A more direct, but substantially more complicated procedure involves the solution 
of transport equation for the PDF’s of the subgrid scalars rather than assuming their 
form. This approach, like its counterpart in turbulence modeling, has the advantage 
that the effect of scalar-scalar correlations appear in a closed form. However, models are 
needed for the molecular diffusion within the subgrid, and one has to resort to mixing 
models (a subject of current intense research) for providing the closure. The approach 
based on guessed PDF methods is feasible and is within the reach of present day com- 
puters. The approach based on solving the transport equation for the subgrid scale PDF 
might require extensive computational resources (Pope, 1990a). These models must ini- 
tially be developed in a simple flow. A homogeneous flow is an excellent configuration 
for this purpose, and will be discussed in the next section. After the establishment 
of a successful model, it may be utilized for simulating more complicated flows. The 
extent of this utilization is dependent on the available computational resources and on 
the performance of the model in rigorous trials. 

3. DESCRIPTION OF THE PROBLEM 

The subject of our DNS is a two-dimensional homogeneous box flow under the 
influence of a binary chemical reaction of the type A + B — * Products. To impose 
homogeneity, periodic boundary conditions are employed in all directions of the flow, 
identifying the box as a homogeneous member of a turbulent universe. This periodicity 
allows the mapping of all the aerothermochemical variables from the physical domain 
into a Fourier domain, thus allowing the implementation of spectral methods for nu- 
merical simulations. The flow field is assumed to be homogeneous and isotropic, and is 
initialized in a similar manner to that of Passot and Pouquet (1987). This involves the 
superposition of a “random” velocity to a zero mean velocity and also includes random 
initial density, temperature, and Mach number fluctuations. These fluctuations have 
certain energy spectra, and the ratio between the compressible and the incompressible 
kinetic energy can be varied to assess the effects of compressibility. The specification of 
the fluctuating field with a random field is to mimic a “probabilistic” turbulent field in 
the context of “deterministic” DNS. This approach is similar to that followed in previous 
direct simulations of turbulent reacting flows (Riley et al., 1986; Givi and McMurtry, 
1988; McMurtry and Givi, 1989). In a laboratory flow, these fluctuations appear as 
the result of interactions with the surrounding universe. Such interactions do not exist 
in the isolated homogeneous flow considered in our simulations. Therefore, in order to 
introduce the “noise,” which plays a central role in laboratory turbulence, these per- 
turbations are randomly superimposed to initialize the fluctuating field for DNS. The 
generated turbulence field is of decaying nature, i.e. there is no artificial forcing mech- 
anism to feed energy to the small wave numbers. While it is more desirable to have a 
stationary turbulence field to focus only on the behavior of the scalar field (Eswaran 
and Pope, 1988; McMurtry and Givi, 1989), it is not clear how to implement such a 
forcing mechanism in the compressible case without modifying the behavior at low wave 
numbers. 
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The scalar fields are defined to be square waves with the two species out of phase 
and at stoichiometric conditions. The species field is assumed to be dynamically pas- 
sive, in that turbulence influences the consequent transport of the scalar field with the 
neglect of reverse influences. In the non-reacting case, the trace of only one of these 
reactants is considered, whereas in the reacting case, the transport of an appropriate 
Shvab-Zeldovich variable is assumed to portray the reactant’s conversion rate. This is 
possible by assuming an infinitely fast chemical reaction rate, and by assuming that the 
reactants have identical thermochemical properties. In this framework, the effects of 
non-equilibrium chemistry are neglected; the inclusion of such effects are postponed for 
future investigations. 

The computational package employed in simulations is based on the modification of 
a code developed by Erlebacher et al. (1987). This code is based on a spectral collocation 
algorithm with Fourier trial functions. All the variables are spectrally approximated on 
iV 2 collocation points, where N represents the number of collocation points in each of 
the directions. The spatially homogeneous flow evolves in time, and at each time step 
jV 2 defines the sample data size for the purpose of statistical analysis. A third order 
accurate Runge-Kutta finite difference scheme is employed for temporal discretization. 
For a trustworthy simulation, the magnitudes of the Reynolds and Peclet numbers must 
be kept at moderate levels and the size of the time step should be kept small. The code 
is capable of simulating both two- and three- dimensional flows and we have performed 
both such simulations. In the presentation of our results in the next section, however, 
we limit our discussion to that of a two-dimensional flow. 

4. PRESENTATION OF RESULTS 

Computations are performed on a domain with a normalized dimension of 2 tt in each 
of the directions of the flow. With the available computational resources, a resolution 
of 256 x 256 collocation grid is attainable. With this resolution, the magnitude of the 
Taylor microscale Reynolds numbers that could be simulated is in the range, Re\ « 
20 30. Simulations are performed with a wide spectrum of compressibility levels; 

here we only report the results obtained by use of two extreme cases: one with a low 
compressibility level, or pseudo incompressible, and the other with a relatively high level. 
In the former, the initial value of the normalized density rms is very small, i.e. p rm » = 
yj< p' 2 >/p 0 ss 0, whereas in the latter, this value is of order unity. The compressibility 
level was monitored by adjusting the initial values of the following parameters (Passot 
and Pouquet, 1987): (1) the rms of the Mach number, M t = 'J < A /' 2 >, and (2) 
the ratio of the compressible energy to that of the total kinetic energy, x- Below, we 
limit the discussion to the results obtained for two cases. In what follows, we refer 
to pseudo incompressible simulations in which the following values are used initially: 
M t = 0.2, x — 0.01, and to compressible simulations in which at the initial time 
M t = 0.6, x = 0.2. 

The compressibility effects can be manifested by both flow visualization and by 
considering the global behavior in an integral sense. In the former, the contour plots 
of the relevant variables constructed from the DNS results show the qualitative behav- 
ior, whereas in the latter, the ensemble averages of the simulated results portray the 
quantitative response. To demonstrate this, in Figs. 1 and 2 we present the contour 
plots of the density for the pseudo incompressible case and for the compressible case, 
respectively. A prominent difference between the two cases is the formation of steep 
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gradients in the high compressible case which are not observed in the pseudo incom- 
pressible simulations. The regions of high gradients are referred to as shocklets, and 
consistent with the previous simulations of Passot and Pouquet (1987), appear when 
the initial levels of density and Mach number fluctuations are high. The effects of in- 
creased compressibility can be quantitatively demonstrated by examining the temporal 
variations of the maxima and minima of the divergence of the velocity, the fluctuating 
Mach number and the ratio of the compressible to the total kinetic energy. These are 
presented in Figs. 3-6. In Figs. 3 and 4 the temporal variations of the minimum and 
maximum values of the divergence of the instantaneous velocity (V • Y) are presented 
for the pseudo incompressible and for the compressible cases, respectively. Note that in 
the low compressible flow, the minimum and maximum values are approximately mirror 
image of each other (with respect to V • Y = 0). In the high compressible case, however, 
the maximum and minimum divergence curves are asymmetric. The minimum value of 
divergence occurs at the normalized time of t « 0.7. This time corresponds to the onset 
of formation of the shocklets, and indicates severe compression within the flow. It must 
be mentioned that the capture of these shocklets by means of global numerical methods 
(without any artificial dissipation) is rather difficult. With 256 x 256 collocation points, 
this was the strongest shock that we were able to capture, and a lower resolution would 
result in a significant oscillation in the magnitude of the velocity divergence after the 
appearance of the shock. One must be careful not to confuse these numerical oscillations 
with compression and expansion waves. 

The effects of compressibility are further quantified in Figs. 5-6. For the pseudo 
incompressible flow, Fig. 5, the magnitude of M t decreases monotonically, whereas in the 
compressible case, Fig. 6, the decrease in M t is interrupted by plateaux. The locations 
of these plateaux coincide with those corresponding to a rise in the compressible kinetic 
energy and the local minima of the velocity divergence. Therefore, these locations 
correspond to the times at which shocklets can be formed in the flow. In the pseudo 
incompressible flow, the magnitude of x remains fairly constant, indicating that the 
initial low level of compressibility remains low at all times. 

With the development of the flow, the species field would consequently evolve. To 
visualize the flow, the contour plots of species A are presented in Fig. 7. Parts (a) 
and (b) of this figure correspond to the zero and the infinitely fast reaction rate cases, 
respectively. This figure exhibits the effects of random motion on the distortion of the 
scalar field and the mixing of the two initially segregated reactants (the contours form 
parallel lines at t = 0). The effects of chemical reaction are to increase the steepness of 
the scalar gradients and to reduce the instantaneous values of the reactants, as indicated 
by a comparison between parts (a) and (b) of the figure. The quantitative behavior of 
the scalar development is depicted in a statistical sense by means of examining the 
evolution of the PDF’s of the conserved Shvab-Zeldovich variable J . This variable is 
normalized and defined within the region [0, 1]. Correspondingly, its PDF, V\{J) is 
always bounded in this region. The temporal variation of Vi(J) is presented in Fig. 8. 
It is shown in this figure that at the initial time, the PDF is approximately composed 
of two delta functions at J = 0,1, indicating the two initially segregated reactants A 
and B. At later times, it evolves through an inverse-like diffusion in the composition 
space. The heights of the delta functions decrease and the PDF is redistributed at 
other J values in the range [0, 1]. At even later times, the PDF becomes concentrated 
around the mean value. Proceeding further in time results in a sharper peak at this 
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mean concentration, and the PDF can be approximated by a Gaussian distribution. 
This Gaussian type behavior has been observed in previous simulations of Givi and 
McMurtry (1988), and Eswaran and Pope (1988), and also has been corroborated in a 
number of experimental investigations. 

An interesting character of these PDF’s is that throughout their evolution, the 
simulated results compare remarkably well with that of a Beta density. This is also 
shown in Fig. 8 by a comparison between the Beta density and the DNS generated 
PDF’s. The Beta density is parameterized with the same first two moments as those 
of the DNS. Therefore, in all the figures, the results are presented with respect to a 
time scale (t*) proportional to the decay of the variance of the scalar J7. Higher order 
moments of the DNS data for the variable J also show good agreements in comparison 
with those predicted by the Beta density. This is demonstrated in Figs. 9 and 10, 
where the normalized kurtosis (^ 4 ) and superskewness (/i 6 ) are presented of the random 
variable J and the reactant A, respectively. At time zero, these moments are close to 
unity and monotonically increase as mixing proceeds. For the Shvab-Zeldovich variable 
J, the magnitudes of the kurtosis and superskewness resulting from the Beta density 
asymptotically approach the limiting values of 3 and 15, respectively. These correspond 
to the normalized fourth and sixth moments of a Gaussian distribution as the variance 
of J tends to zero (as t* -► oo). The DNS generated results are very close to those of the 
Beta distribution throughout the simulations. However, the limiting value for variance 
of J approaching zero can not be obtained in the simulations due to obvious numerical 
difficulties. In the reacting case, the moments of the scalar A are consistently higher 
than those of the conserved scalar, but portray similar trends in both Beta density and 
DNS generated results. 

The trends shown above are also observed in the compressible simulations. The 
profiles of the PDF’s and those of the higher order moments and their comparisons with 
the corresponding parameterized Beta density are presented in Figs. 11-13. Again, the 
comparison is remarkably good. The main difference between these results and those of 
pseudo incompressible simulations is the time scale of the decay of the variance by which 
the PDF is evolved.f This time scale can not be incorporated into the PDF description 
in the format employed here. The single-point nature of these PDF’s do not allow 
for any information about the length scales (or any other scales requiring two-point 
statistics) of turbulence. For a systematic inclusion of the length scale into the PDF 
and quantitative description of the differences between the results in the two cases, one 
is required to consider the evolution of the PDF s at two points (at least). Employing 
the results of DNS to evaluate (or to generate) models for more than one-point-level 
closures (O’Brien, 1985; Eswaran and O’Brien, 1989) has proven to be a challenging 
task, and is the subject of current research (Jiang and Givi, 1991). 

The results of these simulations indicate that the approximation of a Gaussian PDF 
for the final stages of mixing of a conserved scalar is well justified. This corroborates the 
results of laboratory experiments (e.g. Miyawaki et al., 1974; Tavoularis and Corrsin, 
1981) and those of other numerical simulations (e.g. Eswaran and Pope,. 1988; Givi 
and McMurtry, 1988). The numerical experiments performed here, however, are similar 
to most of the laboratory investigations in that they include the results of only one 

| This is not illustrated in the figures clearly, since the time is scaled by the rate of 
variance decay. 
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experimental realization. Future numerical experiments using different initializations 
and/or with better numerical resolution and including three dimensional effects would 
be logical extensions of this work. 

In addition to the Beta distribution approximations, the DNS results also compare 
very well with the mapping closure recently introduced by Kraichnan (see Chen et al, 
1989 and Kraichnan, 1990). This closure has the property of allowing the relaxation of 
the PDF of a conserved scalar property to an asymptotic Gaussian distribution. Pope 
(1990b) has utilized this closure for the prediction of PDF evolution in an isotropic 
homogeneous turbulent flow with an initial double delta distribution (similar to the 
condition imposed here). He showed that in the context of a one point PDF, the results 
obtained by this closure compare favorably with the DNS results. However, again there 
is no length scale information built into the model, and only the evolution of the PDF 
without any information regarding the time scale of evolution can be predicted. Madnia 
and Givi (1991) have extended this model for predicting the decay rate of a reacting 
scalar in homogeneous turbulence. The generated results obtained by applying this 
model for the prediction of the decay rate of a reacting scalar for the cases considered 
in DNS are shown in Figs. 14 and 15. In these figures, the DNS generated results and 
those predicted by a Beta density Eire also shown. The two PDF’s (by the mapping 
closure, and the Beta density) are constructed in such a way to yield the same first two 
normalized moments of the Shvab-Zeldovich variable as those of DNS. In this context, 
the predicted results of the decay by both models show good agreements. The agreement 
for the mapping closure model is rather surprising, considering that this closure was 
developed for an isotropic, fully three-dimensional turbulent flow. This is probably 
due to matching of the normalized first two moments. A better test of the model 
would be its utilization at two-point level, and its comparison with DNS generated 
results. Also, the agreement between the mapping closure results and those of DNS 
worsens as the compressibility level is increased. This again is understandable in that 
Kraichnan’s model was developed for a purely incompressible flow without including 
any compressibility effects. 

With the construction of the DNS database, the results are used to construct the 
PDF’s of the scalar quantities within the subgrid. Since the major riddle in PDF 
modeling is associated with the closure of the diffusion term and not the chemical 
reaction term, we limit the analysis to that of the generated PDF’s of the conserved 
Shvab-Zeldovich variable. The construction of the PDF’s is implemented by placing a 
coarse (n 2 ) mesh over the physical space occupied by the fine mesh (N 2 ). This mesh can 
be considered as the one to be potentially used in LBS. In this way, it is assumed that 
the LES predicts the averaged results at the centers of the n 2 mesh, and the fluctuations 
within the coarse grids are to be modelled by an appropriate PDF closure. The actual 
LES in this case would correspond effectively to simulations with a homogeneous box 
filter in which the values of the filtered means Eire constant. This is somewhat analogous 
to the procedure followed by Schumann (1989). The main difference is the mechanism 
by which the fluctuations are to be considered. Schumann (1989) simply neglected 
such effects, but with consideration of the PDF, it may be possible to account for such 
effects, albeit statistically. Within each cell of the coEirse mesh, we have the values of 
the scalar quantities at N+ = ( N/n ) 2 equally spaced points. This means that there 
is an ensemble of N# sample points at each location to construct the PDF s. This 
construction is implemented by statistical sampling of N $ number of data points. The 
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domain of the scalar property <j> 6 ^mai] is divided into a number of bins of size 

A<£, and the number of scalar variables within each of these bins are counted. The PDF 
of the variable <f>, denoted by V(0 is calculated in the interval ( < <f> < ( + df, d£ = d<f>, 
from the definition of the probability: 

■p(£)<f£ = Probability (£ < <f> 5: £ + d£) (-') 

which translates into: 

V(£)d£ = Probability Density 

Number of elements in the inter val £ < <l> < £ + d( (12) 

= N * 

The value of N+ is a measure of the width of the filter, in that it determines the scale 
at which the fluctuations of the scalar field are not accounted for deterministically, 
but are included in a probabilistic manner by the PDF’s. Consequently, the shapes 
of these distributions are dependent on the magnitude of the sample size N#. If this 
value is too small, then there will be a large scatter in the data. If it is too large, then 
the PDF’s correspond to those appropriate for Reynolds averaging and not LES. In 
this primeval implementation, we have not yet performed an actual LES (using a PDF 
subgrid closure). Rather, we have focused only on constructing the shapes of the PDF’s 
from the DNS results. These PDF’s can again be quantified by the magnitudes of their 
moments to provide a reasonable a priori assessment of the closure for future actual 

LES. 

With the 256 2 available total data, there is a limit on the magnitude of the ensemble 
data for an a priori assessment. Our experience indicated that an ensemble of about 
32 2 was the lowest acceptable limit for a meaningful statistical analysis. Therefore, the 
domain was broken into 8x8 squares, within each the PDF’s were constructed. There 
was also some analysis with 4x4 squares for better statistics. However, since this is 
very close to the Reynolds averaging limit, the generated results will not be considered 
in the forthcoming discussions. 

The PDF’s constructed from this ensemble showed that as the mixing proceeds, all 
the PDF’s tend to have an asymptotic Gaussian distribution. This is to be expected, 
since a subset constructed from a large random set with a probability distribution tends 
to portray the same statistical behavior. However, the results showed that for all the 
cases considered, the PDF’s tend to have an evolution similar to that of the Beta density. 
This is demonstrated in a sample plot of the distributions in Fig. 16 for both pseudo 
incompressible and compressible simulations. In this figure, the PDF’s within one of 
the coarse grids are compared with that of a Beta density with the same two first 
moments as those within the subgrid. The comparison is fairly good, considering the 
size of sampled data. The overall deviation from the Beta density and the Gaussian 
distribution in all the subgrids can be quantified by measuring the norms of the deviation 
of higher moments. The L 2 norms of the fourth moments for the pseudo incompressible 
and the compressible simulations are shown in Figs. 17 and 18, respectively: Note that 
at the initial time, the deviation from a Beta distribution is very small and increases 
with progress in time, but not substantially. At final times, the DNS data are closer 
to that of a Gaussian distribution, but the magnitudes of the second moments are not 
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small enough for the Beta density to approach a Gaussian distribution. Similar overall 
trends are also observed for the compressible simulations. 

With the results of the analysis, albeit primitive, it is speculated that the approach 
based on PDF parameterization from its first two moments may prove serviceable for 
LES, at least for a conserved scalar. However, the approach requires the accurate 
input of the first two moments of the subgrid which must be provided by modelled 
LES transport equations. In the absence of a better alternative, the Beta density may 
prove useful, and the generated results can be used to estimate the maximum reaction 
conversion rate in the infinitely fast chemistry limit. However, the generalization to 
include finite rate kinetics may not be very straightforward, since the first two moments 
(including the covariance of the scalar fluctuations) must be known (or modelled) a 
priori. By the same token, the extension to more complex kinetics with the inclusion 
of non-equilibrium effects, which are very important for practical applications, remains 
to be a challenging task. 

An improvement of the procedure obviously involves the solution of a transport 
equation for the PDF (rather than its parameterization with its moments). However, 
the computational burden may be prohibitive due to the increased dimensionality of 
the subgrid PDF transport equations. An estimate of the computational requirements 
indicates that the cost associated with the implementation of LES procedure involving 
the stochastic solution of the PDF transport equation may be of the same order as 
that of DNS on the fine grids, unless the ratio of the fine to coarse grid is large. In 
this regard, the tradeoff between LES and DNS depends on other factors, such as the 
physical complexity and computational resources. 

5. CONCLUDING REMARKS 

A spectral collocation algorithm has been employed to study the mechanism of 
mixing and reaction in a non-premixed, decaying two-dimensional homogeneous tur- 
bulent flow. The evolution of the species field in a one-step stoichiometric reaction of 
the type A + B — * Products was the subject of the investigation. Calculations were 
performed for zero-rate and infinitely fast rate chemistry in both pseudo incompressible 
and relatively high compressible flows. The effects of compressibility are delineated by 
the formation of thin shocklets and can be quantified by a global examination of the 
data. Mixing and the reaction conversion rate (in the reacting case) are characterized by 
examining the evolution of a Shvab-Zeldovich conserved scalar quantity extracted from 
the DNS. The results show that the PDF of this scalar variable evolves from an initial 
double-delta function distribution to an asymptotic shape which can be approximated 
by a Gaussian distribution. During this evolution, a Beta density qualitatively describes 
the DNS generated PDF’s in both pseudo incompressible and compressible simulations. 
This is quantitatively demonstrated by a good agreement between the values of the 
higher order moments of the PDF’s calculated from the DNS data and those predicted 
by the Beta approximation. The generated results also compare very well- with those 
predicted with the mapping closure of Kraichnan. The decay of the mean of the reacting 
scalar determined by DNS compares very well with that predicted by this model and 
that obtained via the use of the Beta density. This trend is observed in both pseudo 
incompressible and compressible simulations. The main difference is the time scale of 
the decay of the variance of the PDF’s. 
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The Beta density seems to also predict the behavior of the fluctuating field within 
the subgrid for a conserved scalar quantity. Therefore, it may provide a reasonable 
means of parameterizing the subgrid fluctuations in a stochastic sense. The procedure, 
however, requires the input of the first two moments of the variable within the subgrid, 
and these must be provided by modelled LES transport equations. Also, the gener- 
alization to include finite rate kinetics may not be very straightforward, since these 
moments must be modelled a priori. An improvement of the procedure involving the 
solution of a transport equation for the PDF is possible, but may not be practical due 
to the increased dimensionality of the modelled subgrid PDF transport equations. From 
the computational point of view, it is recommended to assess the applicability of such a 
procedure in (yet) simpler flows, before its implementation in simulating more complex 
flows. 

In comparing the DNS results with those of predictions, the first two moments of the 
Shvab-Zeldovich variable must be matched. This is because single-point PDF’s do not 
include information about the frequency scales of turbulence. Therefore, the evolution 
of the length scale, or any other parameter requiring two-point information, must be 
introduced in an ad hoc manner (here, such information is provided via DNS). Future 
validation studies of the PDF’s by DNS need to consider the evolution of multi- (at least 
two-) point statistics. O’Brien (1985) has already employed a two-point formulation in 
the statistical description of homogeneous turbulent flows. The preliminary results 
obtained by such formulations are in good agreement with experimental data (Eswaran 
and O’Brien, 1989). However, a comparison between the predicted results and the DNS 
data, similar to the ones reported here for single point PDF’s is recommended. 

Future work should also deal with finite rate chemistry with the inclusion of non- 
equilibrium effects, and also on examining the extension of PDF methods in dealing 
with the subgrid closures in LES of such flows. 
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Figure 1. Plot of density contours for the 
pseudo incompressible case, t = 6.142. 



Figure 2. Plot of density contours for the 
compressible case, t = 0.CS6. 





Figure 4. Temporal variation of the minimum and maximum values of 
the velocity divergence for the compressible case. 








Figure 7. Plot of species A concentration contours at 
t* = 1.0715. (a) non-reacting case; (b) reacting case. 



Figure 8. Temporal evolution of Vi(J) for the pseudo incompressible 
case. DNS data (dotted line), Beta density (solid line). 

(a) t* = 0.035, (b) t * = 0.549, (c) t m = 0.800, 

(d) r = 1.0715. 















Figure 9. Temporal variation of the kurtosis, ^ 4 , for the 
pseudo incompressible case, (a) variable J , 

(b) Species A in reacting case. 







Figure 10. Temporal variation of the superskewness, /i6, for 
the pseudo incompressible case, (a) variable J , 

(b) Species A in reacting case. 
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Figure 14. Mean concentration of the reacting species A 
vs. time for the pseudo incompressible case. 



t* 


Figure 15. Mean concentration of the reacting species A 
vs. time for the compressible case. 




Figure 1C. Sample PDF’s of the variable J within the subgrid, 
(a) pseudo incompressible case, t * = 0.0355, (b) pseudo incompressible case, t 
(c) compressible case, t * = 0.101, (d) compressible case, t* — 1.507. 
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Abstract. Closed form analytical expressions are obtained for predicting the limiting rate of mean 
reactant conversion in homogeneous turbulent flows under the influence of a binary reaction 
of the type F + rO -> (1 + r) Product . These relations are obtained by means of a single- 
point Probability Density Function (PDF) method based on the Amplitude Mapping Closure 
(Kraichnan, 1989; Chen et ai, 1989; Pope, 1991). It is demonstrated that with this model, the 
maximum rate of the mean reactants’ decay can be conveniently expressed in terms of definite 
integrals of the parabolic cylinder functions. For the cases with complete initial segregation, it is 
shown that the results agree very closely with those predicted by employing a beta density of the 
first kind for an appropriately defined Shvab-Zeldovich scalar variable. With this assumption, the 
final results can also be expressed in terms of closed form analytical expressions which are based 
on the incomplete beta functions. With both models, the dependence of the results on the 
stoichiometric coefficient and the equivalence ratio can be expressed in an explicit manner. For a 
stoichiometric mixture the analytical results simplify significantly. In the mapping closure these 
results are expressed in terms of simple trigonometric functions. For the beta density model they 
are in the form of gamma functions. In all the cases considered, the results are shown to agree 
well with data generated by Direct Numerical Simulations (DNS). Due to the simplicity of these 
expressions and because of nice mathematical features of the parabolic cylinder and the 
incomplete beta functions, these models are recommended for estimating the limiting rate of mean 
reactant conversion in homogeneous reacting flows. These results also provide a valuable tool in 
assessing the extent of validity of turbulence closures for the modeling of unpremixed reacting 
flows. Some discussions are provided on the extension of the models for treating more complic- 
ated reacting systems, including realistic kinetics schemes and multiscalar mixing with finite rate 
chemical reactions in more complex configurations. 

Nomenclature 

a , b y c : some constant. Da: the Damkohler number. 

B: the beta function. F: fuel. 

the parabolic cylinder function. the parameter in the mapping closure. 


1 This work was sponsored by the NASA Langley Research Center under Grant NAG- M 122, and by the Office of Naval 
Research under Grant N00014-90-J-4013. Computational resources were provided by NCSA at the University of Illinois. 

79 


PP€CEDfN€ PAGE B'.ANK NOT FILMED 



80 


C\K. Madnia, S H Frankel, and P Givi 


if' 

H: 

O: 


A. ft 


half the inverse of the normalized variance 


the PDF 

of the Shvab Zeldovich variable. 

r 

the stoichiometric coefficient. 

the Heaviside function. 

t 

the physical time. 

the incomplete beta function. 

W 

the area weight of the reactant. 

the Shvab Zeldovich variable. 

y 

the dummy variable of integration. 

the oxidizer. 

2 

the unmixedness ratio. 

Greek 

Letters 

, . . . : parameters of the beta density. 

<Po 

the composition space for a 

y: the equivalence ratio. 


reference field. 

T: the gamma function. 

X 

the normalized time. 

<5: the delta function. 

X 

the mapping function. 

£: the composition space for the PDF. 




Gaussian 


Subscripts 

0: time zero (inlet of plug (low reactor). st: stoichiometric. 

G : Gaussian. 


Other symbols 

< ): the probability average. the fluctuation for the ensemble mean. 


1. Introduction 

For the past 40 years, since the early work of Hawthorne et al (1949), estimation of the mean 
reactant conversion rate has been the subject of wide investigations in mathematical modeling of 
turbulent reacting flows. In unpremixed reacting systems, including diffusion flames, there are two 
factors by which this rate is influenced: (1) the speed at which the reactants are brought into the 
reaction zone, and (2) the rate at which they are converted to products through chemical reactions. 
The relative importance of the two mechanisms is characterized by magnitude of the Damkohler 
number (Da), which is the ratio of the characteristic frequency of the chemical reaction to that of the 
hydrodynamics. The role of the Damkohler number in the characterization of reacting flows is very 
important (Williams, 1985). In the limit of Da — ► oo, the rate of reactant consumption is governed by 
the hydrodynamics, i.e., the reaction is “mixing controlled” and is determined by the speed at which 
the reactants are brought into an infinitely thin reaction zone (Bilger, 1980). Obviously, with the 
assumption of an infinitely fast chemistry, it is not possible to account for many interesting issues 
associated with nonequilibrium effects in unpremixed flames (Libby and Williams, 1980). However, as 
indicated in the original pioneering work of Toor (1962), and later by O’Brien (1971) and Bilger 
(1980), it is very important to have a prior estimate of the “limiting” rate of mean reactant conversion 
in practical modeling of combustion systems. In this limit the problem reduces to the simpler problem 
of “mixing,” in which its analysis is much simpler (Toor, 1962, 1975). 

Development of an appropriate turbulence model which can predict the mean rate of reactant 
conversion has been the subject of extensive investigations (for reviews see Toor, 1975; Brodkey, 1981; 
Libby and Williams, 1980; Williams, 1985). Amongst the theoretical tools developed, it is now firmly 
accepted that the approach based on the single-point Probability Density Function (PDF) of the 
scalar quantities is particularly useful, and this approach has been very popular for modeling the 
reactant conversion in a variety of turbulent reacting flow systems (Kollmann, 1990; Pope, 1985, 1990, 
1991). The advantage of PDF methods is due to their inherent capability to include all the single- 
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point statistical information pertaining to the scalar field. Therefore, once the PDF (or the joint PDF) 
of the scalar variables is determined, all the relevant one-point statistics of the Held are available 
without need for additional closures. 

The most logical and systematic means of determining the PDF involves the solution of an 
appropriate transport equation governing its evolution. In this equation, due to the nature of the 
formulation, the effects of chemical reaction appear in a closed form (Pope, 1976), regardless of its 
degree of complexity. However, the influences of molecular action cannot be fully described, and can 
be treated only by means of employing an appropriate closure. As noted by Pope (1991), in most 
previous applications this problem has been circumvented through the use of the Coalescence/ 
Dispersion (C/D) models. Examples of such models are the early C/D prototype of Curl (1963), the 
Linear Mean Square Estimation (LMSE) theory (O’Brien, 1980), the closure of Janicka et ai (1979), 
among others (Pope, 1982, 1985; Kosaly and Givi, 1987; Givi, 1989; Dutta and Tarbell, 1989). 
Despite their advantageous characteristics, the shortcomings associated with the C/D closures in the 
probabilistic description of scalar transport are well recognized. Namely, none of the aforementioned 
models predict an asymptotic Gaussian distribution for the PDF of a conserved scalar variable in 
homogeneous turbulence (Pope, 1982, 1991; Kosaly and Givi, 1987); and those which are capable of 
doing so (e.g., Pope, 1982), do not predict the initial stages of mixing correctly (Kosaly, 1986; Kosaly 
and Givi, 1987). 

Recent development of the Amplitude Mapping Closure by Kraichnan and coworkers (Kraichnan, 
1989; Chen et ai, 1989) (see also Pope, 1991) has provided a promising way of alleviating some of the 
problems associated with the C/D closures. This closure, in essence, provides a means of accounting 
for the transport of the PDF in composition space, and its validity and physical applicability 
have been evidenced in a number of comparisons against data generated by means of both direct 
numerical simulations (Pope, 1990, 1991; Gao, 1991a, b, Madnia et ai , 1991b; Jiang et ai , 1992) and 
laboratory experiments (Frankel et ai , 1992a). These results suggest that, at least in the setting of an 
isotropic turbulent flow, this closure has some superior features over all the previous C/D-type 
models. 

Based on this demonstrated superiority, our objective here is to examine further the properties of 
this closure and to assess its capabilities for applications in modeling of unpremixed turbulent 
reacting flow systems. In particular, it is intended to provide a reasonably simple recipe that can be 
used in conjunction with this closure for predicting the limiting rate of mean reactant conversion. 
However, since this is the first study of this type, and due to mathematical complexities (that soon 
become apparent), we have made some simplifying assumptions which are indicated here. Firstly, we 
consider an idealized irreversible binary reaction of the type F + rO-^(l + r) Product with initially 
segregated reactants (F and 0 ). In accordance with the discussions above, only the maximum rate 
of mean reactant conversion is considered. Secondly, the turbulence field is assumed statistically 
homogeneous. Thirdly, all the chemical species are assumed to have identical and constant thermo- 
dynamic properties. Finally, the flow field is assumed isothermal in which the dynamic role of the 
chemical reaction on the hydrodynamic field is ignored. 

With all these assumptions, the reacting system considered is obviously an idealized prototype of 
conventional combustion systems. However, it does provide a good model for dilute reacting systems 
in typical mixing controlled plug flow reactors (Toor, 1962, 1975; Bilger, 1980; Hill, 1976; Brodkey, 
1981). Moreover, because of the mathematical complexities, even in this simple case, it is deemed 
necessary to analyze this simplified system before considering more complex scenarios. Nevertheless, 
the model is capable of accounting for arbitrary values of the stoichiometric coefficients and for any 
equivalence ratio. This allows the capture of many interesting features, as will be demonstrated. 

For the idealized case of initially segregated reactants, the initial marginal PDFs of their 
concentrations are composed of “delta functions.” Therefore, it is speculated that the approach based 
on an assumed probability distribution may also provide a reasonably good closure. Therefore, in 
addition to the mapping closure, a member of the family of Pearson frequencies is also considered. 
The results obtained by this frequency are compared with those of the mapping closure, 'and are also 
assessed against data generated by means of DNS. 

In the next section the problem under consideration is outlined along with the mathematical basis 
by which the single-point PDF methods are used. In Section 2.1 the salient features of the mapping 
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closure at the single-point level are analyzed including a discussion on the formalities of the closure 
for our purpose. With this closure, a closed form analytical expression is provided for the limiting 
bound of the mean reactant conversion. This limiting rate is also predicted by means of the beta 
density model in Section 2.2. In both these subsections the simplifications for the case of a stoichio- 
metric mixture are made, motivating the use of our final "simple" analytical expressions in practical 
modeling of stoichiometric plug flow reactors. In Section 3.1 the results predicted by both models are 
compared against those generated by DNS of a three-dimensional homogeneous turbulent flow. In 
Section 3.2 some discussions are presented highlighting the implications of these results in turbulence 
modeling. This paper is drawn to a close in Section 4 with some discussions for possible future 
extension of the two models for the statistical description of more complicated chemically reacting 
turbulent flows. 


2. Formulation 

With the assumptions described in the introduction, the statistical behavior of the reacting field in the 
reaction F -f rO -*>(1 + r) Product can be related to the statistics of an appropriate conserved Shvab- 
Zeldovich mixture fraction, / (Bilger, 1980). This mixture fraction can be normalized in such a way to 
yield values of unity in the fuel F stream and zero in the oxidizer 0 stream. For the purpose of 
statistical treatment, we define t ), & { 0 ( t\ and t), respectively, as the marginal PDFs of 

the concentration of F, the concentration of 0, and the Shvab-Zeldovich variable /, For initially 
segregated reactants with no fuel in the oxidizer stream (and vice versa), the initial conditions for the 
marginal PDFs of the concentrations of the two reactants are given by 

*f(1 0) = W f 6(Z - Fj) + W 0 6(a 

^oK. o) = w 0 &(t - o t ) + w F t m o < { < i. {i) 

Here, F t and 0; denote the initial concentrations of the two species in the two feeds, and W F and W Q 
represent the relative weights of the reactants at the initial time (i.e., the area ratios at the inlet of a 
plug flow reactor). With the normalized value of the concentrations equal to unity at the feeds, i.e., 
Fj = 0 { = 1, the stoichiometric value of the Shvab-Zeldovich variable, Xu is determined from the 
parameter r. With the assumption of an infinitely fast chemistry, the marginal PDFs of the reactants 1 
concentrations are related to the frequency of the Shvab-Zeldovich variable (Bilger, 1980; Kosaly and 
Givi, 1987): 

t) = (i - / St )csy/ S , + «i - A), o + s f (tm 

t) = i - a o + s 0 (')<5(£). 

Here, £ > 0 and 

3,(0- 0<*£, 

r r « 

3 o (0 = j 1 -MO. 

The initial condition for the PDF of the Shvab-Zeldovich variable is given by 

&s(Z,0)=W F &tf-l)+W o S(Z). (4) 

Equation (4) implies that <X)(t = 0) = W F . Since X is a conserved variable, its mean value remains 
constant, i.e., <X)(0 — “ Wf* The integration of (2) yields the temporal variation of the 

statistics of the species field at all times, if the PDFs of f are known. As indicated above, in the 
setting of a mixing controlled reaction this PDF provides all the desired statistical properties of the 
reacting field. 


2.1. Amplitude Mapping Closure 

The implementation of the amplitude mapping closure involves a mapping of the random field of 
interest £ to a stationary Gaussian reference field </> 0 , via a transformation £ = *(<Po» 0- Once this 
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relation is established, the PDF of the random variable £, .*(4), is related to that of a Gaussian 
distribution. In homogeneous turbulence, the transport equation for this function satisfies (Chen et ul ., 


1989; Pope, 1991) 


dz 


? 2 x 


+ cVo’ 


(5) 


In this equation, t is a normalized time within which the scalar length scale information is embedded. 
The relations between this time and the physical time, i.e., i(f), cannot be determined in the context of 
single-point PDF description and must be provided by external means (Pope, 1990; Jiang et al , 1992). 
For the case considered here, with the initial PDF of the variable / given by (4), the corresponding 
form of the initial mapping is 

X((po> 0 ) = m<p 0 ~ (/?*)» oc < <p 0 < X), ( 6 ) 

where H is the Heaviside function and is a measure of the initial asymmetry of the initial PDF 
around the ensemble mean of the variable /: 

tp* — s/2 erf _, (l — (7) 

where “erf" denotes the error function. The mapping function is obtained by solving (5) subject to 
initial condition (6). The general solution of this equation has the form (Gao, 1991a) 




X(Vo, *) = 


V 


+ 1 




x(y. 0) exp 




((P 0 e-' - y) 2 (\ + ^ 2 )" 


W 1 


dy. 


where &(z) = ^/exp^r) — 1. Inserting (6) for x{<p 0 , 0) in (8), we have 

X(<Po . t) = i[l + erf(a ^ 0 + &)], 


( 8 ) 

(9) 


where 


a(t) = 


b(z) = 




( 10 ) 


jb' J29 ' 

Finally, the solution for the PDF of / is determined directly from the mapping relation between 
the physical field 4 and the Gaussian reference field q> 0 : 


&/(x(Vo< T X t) = &g(<P o) 


(*lY\ 

Wo/ 


(11) 


Here, 9 G denotes the PDF of a standardized Gaussian distribution, i.e., A(<Po) = ( lis/ln) exp( - <po/2). 
A combination of (11) and (9) yields the final result for the PDF of the Shvab-Zeldovich variable: 


„ , , , , [(<Po e ~ T ~ V*) 2 Vo 

&s(x(<Po> t), t) = ^ exp — _ 1 


( 12 ) 


With a combination of (12) and (2), all the pertinent single-point statistics of the reacting field are 
determined. The most important of these statistics are the ensemble mean values of the reactants’ 
concentrations. These mean values are obtained directly by integrating their respective PDFs. The 
intermediate steps in deriving these relations are not presented but are provided by Frankel (1992). 
Here, only the essential steps are presented. For the mean fuel concentration, <F>, the first part of (2) 
reads 


<F>(t)= P F(QPs(e.T)d( 


-i; 


F(x(<p 0 ' T ). *)&g(Vo) d(p 0 , 


(13) 


where the lower limit of the last integral corresponds to the value of <p 0 at which / is equal to the 
stoichiometric value of the Shvab-Zeldovich variable. Evaluating this limit from (9), equation (13) can 
be analytically integrated. This is possible by representing the error function in ^the form of 
its definition, and performing the resulting definite, double exponential integral. The results after 
extensive algebraic manipulations yield 


<F)(r) 


(1 - 2 /„) 


4(1 - A) 


1 + erf 


b/a + 


- 0 ] 


1 


J2 ) J n^2a(\ - /„) 


«■>(-£ -T-7>~ ««> 
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Similar expressions can be obtained for the mean oxidizer concentration; 


<0>(T) = 


_ erf (^±G] - --i_exp(— 4 - C -- h S 

V s/2 )\ Ky/ '2a/ u \ ^ 2 a 


K'-4 

In these equations, c is related to the stoichiometric coefficient, 

c = — erf 1 (2/ it - 1), 

a 




and 


© + = 


dy expl — a 2 c 2 y 2 + 


A 

Sy 2 




( 15 ) 


(16) 


,tmy + 2S' 


( 17 ) 


__ 2 b c 

y + = + 2acy ----- 
a a 

Here, and are, respectively, the parabolic cylinder functions of order —2 and — 1, 

belonging to the family of degenerate hypergeometric functions (Abramowitz and Stegun, 1972). 

Due to nice properties of the degenerate hypergeometric functions, many of the interesting features 
of (14) and (15) can be depicted. First, simple manipulation of this equation shows that for a 
stoichiometric mixture, = A* both reactants decay at the same rate, i.e., 


<F>(r) _ <0>(r) 

<F>(0) = W F <O>(0) = W Q ’ 


and for a nonstoichiometric initial 
asymptote to 


condition, the limits of the concentration values as t, ■§ -* oo. 


0, 


A > </>, 


Lim <F>(t) = -i 


</> - A 

, l-A ’ 


fo, 


A ^ </>, 


A < </>, 


(18) 


Lim <0 >(t) = 


1 - 


</> 


A > </>, 


These limits are obtained by employing the Taylor series expansion of the relevant functions as 
^ x, and indicate the limiting bound of the concentrations of the unconsumed reactants in both 
fuel-rich and fuel-lean mixtures. While these limiting behaviors are rather trivial from a physical 
standpoint, in a computational procedure it must be made sure that they are satisfied. Because of the 
mathematical properties of the parabolic cylinder functions, these limiting cases can be realized in 
our computational procedure in a relatively easy manner. It would be very difficult to obtain these 
limiting behaviors numerically in an integration procedure within the original unbounded domain. 

At first glance, (14) and (15) may appear somewhat complicated. However, due to nice mathe- 
matical properties of the parabolic cylinder functions (Abramowitz and Stegun, 1972), these 
equations can be integrated rather easily within the finite domain (0 ^ y < 1). For fuel-lean or 
fuel-rich mixtures, the integration can only be done by means of employing numerical methods. 
However, for a stoichiometric mixture, the results simplify further as demonstrated below. 


Stoichiometric Mixture. For practical applications in stoichiometric plug flow reactors, the equations 
simplify considerably. For a stoichiometric mixture, and an initially symmetric PDF around the mean 
value (i.e., <</> = A = i)< both parameters b and c are zero. Under this condition, the first terms on 
the right-hand sides of (14) and (15) drop. Knowing S>_ 2 (0) = i^fO) = 1, the remaining terms yield 


<F) (t) 
<F>(0) 


<0>(t) 

<O>(0) 



dy 


2 y 2 -I- 1/a 2 


1 - 


2 arctan #(r) 


(19) 
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The simplicity of this equation is noteworthy and very pleasing. Because of this simplicity, (19) is 
strongly recommended for engineering predictions of the mean reactant conversion rate in stoichio- 
metric homogeneous flows, such as the plug flow reactors considered in numerous previous investiga- 
tions (Toor, 1962, 1975; Brodkey, 1981; Hill, 1976; O'Brien, 1971; Kosaly and Givi, 1987). 


2.2. The Beta Density Model 


For initially segregated reactants, the initial PDF of the Shvab-Zeldovich variable is composed of 
two delta functions at the extreme limits of the variable. Therefore, it is proposed that the family of 
Pearson (1895) frequencies may provide a reasonable means of estimating this distribution at all times 
(Frankel et «/., 1991; Girimaji, 1991a). The appropriate form of the Pearson distribution is in this case 
the beta density of the first kind. This density has been employed in the statistical description of 
turbulent reacting flows by Rhodes (1975), Jones and Priddin (1978), Lockwood and Moneib (1980), 
Peters (1984), and Janicka and Peters (1982) amongst others (for recent reviews, see Givi, 1989; 
Priddin, 1991). For an initially nonsymmetric PDF, the beta density corresponds to Pearson Type I 
and for the symmetric case to Pearson Type II. The relevance of the latter in modeling of molecular 
mixing from an initial symmetric binary state has been described by Madnia et ai (1991a) and 
Girimaji (1991a). In both cases the PDF of the Shvab-Zeldovich variable is represented by 
(Abramowitz and Stegun, 1972) 

*?>«) = unr -flT ^'" 1(1 - A*' 1 ’ (20) 

B(P i.ft) 

where denotes the beta function, and the parameters ft and ft are dependent on the mean 

and the variance of the random variable # . In applications to the mixing controlled reaction 
considered here, we assume that the PDF of the Shvab-Zeldovich variable always retains a beta 
distribution. Thus all the statistics of the reacting scalar are subsequently determined. The ensemble 
mean values are determined by a combination of (20) and (2). Following the same procedure as that 
described in the section on the mapping closure, after some manipulations the Final results can be 
expressed as (Frankel, 1992) 


<F>(t) = 


/,?■(! -Z./ 2 ' 1 
(ft +P?i.W 


+ 


i ( ft 
i - A Vft + ft 



(i - s,jp lt ft)), 


( 21 ) 


<0>w = 


/» ,rl (i - A/ 2 

(ft + ft)fl(ft,ft) 


+ 



ft 


(ft + ft) A 


•//.(ft.ftX 


( 22 ) 


where ./ denotes the incomplete beta function (Abramowitz and Stegun, 1972). 

Due to nice mathematical properties of the beta function, the final results are cast in terms of its 
integral. In this case, however, the integral can be expressed in terms of the incomplete beta function. 
The mathematical properties of this special function are well known, and the expressions above 
are conveniently amenable to numerical integration (Frankel, 1992). Again, the physical limiting 
conditions discussed before are realized by (21) and (22). That is, in a stoichiometric mixture, both 
reactants decay at the same rate; and in lean or rich mixtures, the same limiting conditions as those in 
(18) are realized. 


Stoichiometric Mixture. Again, in the case of an initially symmetric PDF under stoichiometric 
conditions, the final expressions become simpler. Under this condition, ft = ft, and knowing 
J m (x, x) = }, (21) and (22) reduce to 

<F)(Q _ <0X0 1 r(g) 

<f>( 0) <o>(0) y^r ( ff + i)’ 

where g is half the inverse of the normalized variance of the Shvab-Zeldovich variable, and T denotes 
the gamma function. 

3. Results 


The final forms of (14), (15), (19), (21), (22), and (23) are gratifying since they provide a relatively 
simple and effective means of estimating the maximum rate of mean reactant conversion in homogen- 
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eous reacting flows. As indicated before, the mathematical operations leading to these equations are 
somewhat involved, but the final results can be conveniently expressed in terms of known special 
functions. However, since both models are based on single-point PDF descriptions, these equations 
are not in a complete closed form and are dependent on the parameters g and/or In this context 
this parameter cannot be determined by the model and must therefore be specified by external means 
(Pope, 1991; Jiang et al . , 1992; Frankel et al. t 1992a). This deficiency is not particular to the two 
models considered here, and exists in any single-point statistical description, including all of those 
based on the C/D models. 

The extent of validity of these simple relations can be demonstrated by a comparison between the 
model predictions and the data obtained by means of DNS. The comparison is made here for several 
values of </> and for the purpose of demonstration. In this comparison the magnitudes of the nor- 
malized variance of the Shvab-Zeldovich variable are matched with those of DNS. This implies that, 
for given values of </> and / st , the parameters g , ft, and ft are provided externally from the DNS 

data. With this provision, the model prediction results can be directly assessed against DNS data. 

The DNS procedure is similar to that of previous simulations of this type. For a detailed 
description we refer the reader to Madnia and Givi (1992). The subject of the present DNS is a 
three-dimensional periodic homogeneous box flow under the influence of a binary reaction of the type 
described above. The initial species field is composed of out-of-phase square waves for the two 
reactants F and O. The computational package is based on the modification of a spectral-collocation 
procedure using Fourier basis functions developed by Erlebacher et al. (1990a) (see also Erlebacher et 
al . , 1987, 1990b). The hydrodynamic field is assumed isotropic, and is initialized in a similar manner 
to that of Erlebacher et al. (1990a) and Passot and Pouquet (1987). The turbulent field is of a 
decaying nature in that there is no artificial forcing mechanism to feed energy to low wave numbers. 
The code is capable of simulating flows with different levels of compressibility (Hussaini et al . , 1990). 
Here, only the results obtained for a low compressible case are discussed, since most previous analyses 
of plug flow reactors have dealt primarily with incompressible flows (Toor, 1975; Hill, 1976; Brodkey, 
1981; Leonard and Hill, 1988a, b, 1991). The resolution consists of 96 collocation points in each 
direction. Therefore, at each time step 96 3 is the sample size for statistical analyses. With this 
resolution, simulations with a Reynolds number (based on the Taylor microscale) of Re x ^ 41 are 
attainable. The value of the molecular Schmidt number is set equal to unity. 

3.1. Validations 

The statistical behavior of the scalar field is depicted by examining the evolution of the PDFs of the 
Shvab-Zeldovich variable /. These are shown in Figures 1 and 2 at times close to the initial (tl) and 
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Figure 1. PDF of the Shvab-Zeldovich variable at two times Figure 2. PDF of the Shvab-Zeldovich variable at two times 
(r 2 > /l) for the symmetric case {</> = i). (t2 > 1 1) for the nonsymmetric case (</> = 0.625). 
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Figure 3. Normalized mean concentration of fuel and the 
oxidizer for the symmetric case, and /„ = 0.5. 

the final (f2) states. These figures correspond, respectively, to the cases of an initially symmetric 
(<</> = i) and nonsymmetric = 0.625) PDFs. At the initial time, the PDF is approximately 
composed of two delta functions at J' = 0, 1 indicative of the two initially segregated reactants, F and 
O. At later times, the PDF evolves through an inverse-like diffusion in composition space. The heights 
of the delta functions decrease, and the PDF is redistributed at other / values in the range [0, 1], 
and subsequently becomes centralized around the mean value. Proceeding further in time results in a 
sharper peak at this mean value, and in both cases the PDF can be approximated by a Gaussian 
distribution near the mean scalar value. The trend for the symmetric case is the same as that 
presented in earlier DNS studies (Eswaran and Pope, 1988; Givi and McMurtry, 1988). For the 
nonsymmetric case there are no DNS data in the literature, but the present results verify that the 
asymptotic PDF can still be approximated by a Gaussian distribution near its mean value. 

The PDFs obtained by the mapping closure and those by an assumed beta density are also 
presented in Figures 1 and 2. In these figures the model PDFs are parametrized with the same first 
two moments obtained from DNS. In this parametrization only the normalized magnitude of the 
variance of the models are forced equal to that of the DNS and no attempt was made to account 
for the departure from the “exact” initial double delta distribution in DNS. With this matching, 
nevertheless, the results clearly indicate that the model predictions compare very well with the DNS 
results. Also, both models yield an asymptotic Gaussian-like PDF. 

The temporal variation of the ensemble mean of the reactants’ concentration by the two models 
are compared against those of DNS in Figures 3 and 4. These figures correspond to the two cases of 
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Figure 4. Normalized mean concentration of fuel and the oxidizer for the nonsymmetric case, (a) /„ = 0.4, (b) /„ = 0.2. 
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symmetric and nonsymmetric initial PDFs, respectively. In the symmetric case, under stoichiometric 
conditions, the results are simply obtained from the analytical expressions in (19) and (23). In the 
nonsymmetric case, numerical integration of (14) and (15) and evaluation of the incomplete beta 
function in (21) and (22) are necessary. In both cases the agreement of the models with the DNS data 
is noteworthy. Also, a comparison between parts (a) and (b) of Figure 4 shows that as the magnitude 
of./,, decreases, the rate of consumption of the oxidizer increases. 

The agreements noted above follow from the compatibility of the model PDFs and those of DNS, 
at least for the case of binary mixing considered here. This finding is not new and has been well 
documented in previous works, at least those considering an initially symmetric PDF (Pope, 1990, 
1991, Madnia and Givi, 1992; Madnia et ul., 1991a; Girimaji, 1991a). However, a nice feature of the 
models is the explicit form of the final equations expressing these statistical quantities. Supported by 
this quantitative agreement, it is proposed that, in the absence of a better alternative, the relations 
obtained above be used as an explicit simple means for predicting the maximum rate of mean reactant 
conversion in homogeneous reacting systems. 

Despite the simplicity of these equations and their ease of application, it must be mentioned that 
these equations predict the rate of mean decay of the reactants’ concentration far better than all 
the previous turbulence closures based on the C/D models (see Givi, 1989). This is particularly 
advantageous in that this evaluation can be made by a simple algebraic relation, whereas the C/D 
implementations usually require more expensive numerical simulations (Kosaly and Givi, 1987; 
McMurtry and Givi, 1989). Even for nonunity equivalence ratios, the numerical integration required 
by the two models above is considerably less computationally demanding than those of the C/D 
models. The only input in these models, similar to those in C/D closures, is the variance of the 
Shvab-Zeldovich variable. This is provided here by means of DNS. In an actual implementation, this 
variance can be obtained from experimental data or by means of an appropriate turbulence model 
(Frankel et al., 1992a). The provision of such data is not very difficult since they can be obtained in 
the setting of a nonreacting flow. 


3.2. Applications 

The relations obtained here can be used in determining the extent of validity of other conventional 
closures for predicting the limiting rate of mean reactant conversion in turbulent flows. As an 
example, here we consider the model based on the famous hypothesis of Toor (1962, 1975), which has 
received considerable attention in practical modeling of unpremixed homogeneous reacting systems 
(Bilger, 1980; Brodkey, 1981; Leonard and Hill, 1987, 1988a, b; Kosaly and Givi, 1987- Kosaly 1987- 
Giv! and McMurtry, 1988; McMurtry and Givi, 1989; Givi, 1989). According to this hypothesis in 
an jsothermal homogeneous reacting turbulent flow, the decay of the unmixedness, denoted by 
T - <F O'>(r)/<F'O'>(0), where the prime quantities indicate fluctuation from the ensemble mean 
value, is independent of the magnitude of the Damkohler number. This implies that the normalized 
unmixedness parameter, defined by 




(24) 


is the same under both reacting and nonreacting conditions, i.e., 2(t) = constant = 1 for all values of 
Da. In previous DNS assessments of this hypothesis, it has been shown that for the case of initially 
segregated reactants this model cannot be employed, and the normalized unmixedness ratio depends 
on the nature of mixing and the magnitude of the Damkoher number (Givi and McMurtry, 1988; 
McMurtry and Givi, 1989). In particular, it has been demonstrated that even for Da -* x, while the 
normalized unmixedness is equal to unity at the initial time, its limiting lower bound depends on the 
asymptotic frequency of the Shvab-Zeldovich variable. That is, J"(r = 0) = 1 > _Z(r ) > jZ(r ->x)=C, 
where C is the lower limiting bound. For an asymptotic Gaussian distribution, it can be easily shown 
that, for a mixture under stoichiometric conditions, the lower bound limits asymptotes to the constant 
value C = 2/ zr (Kosaly, 1987; Givi and McMurtry, 1988). 

This deviation from unity can be realized by the two models considered. With the mapping closure, 
under symmetric stoichiometric conditions, from (2), (14), and (15), following the same integration 
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Figure 5. The unmixedness ratio at several values of the 
equivalence ratio. 
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procedure as before, it is shown that 

(FO r y(Da 




2(arctan(l/^)) 2 


(25) 


< F'0'y(Da = 0) n arctan(l/^ v -' / ^ 2 + 2) 

For the beta density model, the corresponding form of the unmixedness ratio for p x = is given by 

(26) 


*v- 


2 gf V(g) V 

* \ng + i)J ' 


In (25) and (26), the subscripts m and P are added to denote the mapping and the beta density model. 
It is easy to show that these equations satisfy the correct limiting conditions for a stoichiometric 
mixture. This is for both the initial time, i.e., the inlet of the reactor, 


Lim = Lim = 1, (27) 

pf-o] [»- 1/2] 

and at large distances from it: 2 

Lim Z m = Lim Z, = C = -. (28) 

The latter limiting condition cannot be realized in any of the previously used C/D models, or by 

means of Toor’s models (McMurtry and Givi, 1989; Givi, 1989). 

The results based on the applications of Toor’s model become less accurate for nonstoichiometric 
mixtures. For equivalence ratios other than unity, with the depletion of one of the reactants, the 
unmixedness parameter approaches zero faster. This is demonstrated by the solution of the mapping 
closure shown in Figure 5 for several values of the equivalence ratio (y). Note that as the magnitude 

of this ratio increases above one, the unmixedness ratio goes to zero more rapidly. For an unity 

equivalence ratio, the correct asymptotic value of 2jn is realized. 


4. Extensions for Modeling of More Complex Reacting Turbulent Flows 

Despite the pleasing features of our simple mathematical expressions, there are several restricting 
assumptions which were necessarily imposed in deriving these equations. Here, we would like to 
address the ramifications associated with these assumptions, and to provide the means of overcoming 
them in future extensions of these models. 

Firstly, due to the assumption of infinitely fast chemistry, only the maximum rate of the mean 
reactant conversion is obtained. While this rate is very useful in describing unpremixed flames, from 
both a theoretical standpoint and for practical applications (Givi, 1989; McMurtry and Givi, 1989; 
Kosaly and Givi, 1987; Toor, 1962, 1975; CTBrien, 1971; Bilger, 1980; Williams, 1985; Dutta and 
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Tarbcll, 1989), the model is not capable of describing some of the important features of the turbulent 
(lames, especially those associated with nonequilibrium effects. The extensions to finite rate chemistry, 
reversible reactions, nonequilibrium flames, and multistep kinetics systems require numerical inte- 
gration of the PDF transport equation. For these cases the problem cannot be mathematically 
reduced to that of keeping track of a single scalar variable (like /), and requires the use of 
multivariate statistical descriptions. For this, the implementation of mapping closure is relatively 
straightforward since it provides a transport equation for the joint PDFs of the scalar variable in a 
multivariable sense (Pope, 1991; Gao and O'Brien, 1991). However, it is not presently clear how to 
devise an efficient computational procedure, typically based on Monte Carlo methods (Pope, 1981), 
for the numerical treatment of these equations. Some work in this regard is currently under way 
(Valino and Gao, 1991; Valino el at., 1991). 

The extension of assumed distributions based on the beta density for treating multiscalars is more 
straightforward but less trivial to justify. The most obvious means is to implement the multivariate 
form of the Pearson distributions. The direct analogue of the beta density is the Dirichlet frequency 
(Johnson, 1987, Johnson and Kotz, 1972; Wilks, 1962; Narumi, 1923). For a mixture composed of 
•V + 1 species, the joint PDF of N concentrations (ip t , i/» 2 , .... <j/ y ) is described in terms of an 
,\ r -variate density of the form 


■Ay 


_ r(/?i + P 2 + "' + Av+i ) 

r(A m ),..., n/? v+l ) 




- l/b - ~ 


(29) 


subject to the physical constraint 

I = !• (30) 

The application of this density in modeling of multispecies reactions has been nicely discussed by 
Girimaji (1991a). Due to the mathematical properties of the gamma function, this density is pleasing 
from a mathematical viewpoint and most statistical cross correlations of the random variables 
(•/'i, •• •) can be conveniently obtained by means of simple analytical relations (Frankel, 1992). 

Some points in this regard have been made by Girimaji (1991b). However, the use of the Dirichlet 
frequency cannot be justified for modeling of unpremixed. reacting flow in a general sense (Frankel, 
1992). In fact, there is no way of implementing this density directly for modeling of nonequilibrium 
fiames, involving strong temperature variations. This is simply due to the additivity constraints of this 
density requiring the unity sum of the normalized random variables (30). 

Secondly, the mathematical derivations presented here are only valid for initially segregated 
reactants. In both models the complete segregation facilitates significant simplifications of the final 
equations. This assumption is compatible with that made in the majority of previous works on 
unpremixed reacting flows (Toor, 1975; Brodkey, 1981; Bilger, 1980). For other initial conditions, 
e 'g', P artia * premixing of the reactants, or non-delta-like distributions, numerical integration of 
the PDF transport equation is required. Again, an appropriately devised numerical procedure can 
accommodate such conditions. However, the use of a beta density (or any other assumed distribu- 
tions) cannot be justified for other complex initial conditions. 

Thirdly, the final mathematical expressions are only valid in the setting of a homogeneous flow. 
Extension to inhomogeneous flow predictions is also straightforward, but requires numerical inte- 
gration of the modeled equations. Both models can be directly implemented into appropriately 
devised numerical procedures. The mapping closure can be invoked in the mixing step of a fractional 
stepping procedure, similar to that of typical Monte Carlo procedures (Pope, 1981). The beta density 
requires modeled transport equations for the low-order statistics of the reacting field. These equations 
include the required information pertaining to the spatial inhomogeneity of the flow through the 
parameters /?,, .... With this information, all the higher-order statistics of the reacting field can be 

provided by simple analytical means (Girimaji, 1991b; Frankel et al., 1992b). 

Finally, in the context of the single-point PDF formulation presented, there- is no information 
pertaining to the evolution of the relevant turbulent length scales. The final expressions can only be 
presented in terms of other physical parameters (here, through the variance of the Shvab-Zeldovich 
variable). In the context considered, this parameter has been provided by the DNS data. In a practical 
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application, this must be provided by external means (turbulence models, experimental data, etc ). 
Jiang et al. (1992) and Frankel et al (1992a) provide further discussions related to this issue. 


5. Concluding Remarks 

It is demonstrated that the mapping closure of Kraichnan (Chen et al . , 1989; Pope, 1991) yields 
closed-form analytical expressions for predicting the limiting bound of the mean reactant conversion 
rate in simple chemistry of the type F + rO -►(! + r) Product in homogeneous, isothermal turbulent 
flows. It is also shown that, for the case of complete initial segregation, the scalar PDFs generated by 
this closure can be well approximated by a beta density. This density also provides closed-form 
analytical expressions for the limiting rate of mean reactant conversion. A nice feature of the 
mathematical results generated by the two models is their capability of revealing the influence of the 
stoichiometric coefficient and the equivalence ratio. In both cases the mathematical expressions 
simplify significantly for a stoichiometric mixture. The prediction results via both models compare 
favorably with data generated by DNS. This agreement follows from the compatibility of the models’ 
PDFs with those of DNS. The simple final results generated here are superior to those of previous 
closures based on typical C/D models, and those based on Toor’s hypothesis. 

These closed-form relations are furnished with the imposition of several restrictive assumptions. 
The ramifications associated with these assumptions are discussed, and some suggestions for future 
extensions are provided. Despite these assumptions, it is very encouraging to have physically plausible 
algebraic relations for the direct estimate of the mean reactant conversion rate in homogeneous 
turbulent flows, typical of those in plug flow reactors. 
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Results are presented of direct numerical simulations (DNS) of a spatially developing shear flow 
under the influence of an infinitely fast chemical reactions of the type A + B-* Products. The 
simulation results are used to construct the compositional structure of the scalar field in a statistical 
manner. The results of this statistical analysis indicate that the use of a Beta density for the probability 
density function (PDF) of an appropriate Shvab-Zeldovich mixture fraction provides a very good 
estimate of the limiting bounds of the reactant conversion rate within the shear layer. This provides a 
strong justification for the implementation of this density in practical modeling of non-homogeneous 
turbulent reacting flows. However, the validity of the model cannot be generalized for predictions of 
higher order statistical quantities. A closed form analytical expression Is presented for predicting the 
maximum rate of reactant conversion in non-homogeneous reacting turbulence. 

KEYWORDS Direct numerical simulations Turbulent shear flows Beta density 
Probability density function Unmixedness. 


INTRODUCTION 

In our previous work (Madnia et al., 1991a) we have demonstrated that the Beta 
density provides a very good means of approximating the probability distribution 
of a conserved scalar quantity in homogeneous turbulent flows. The versatility of 
this density enables it to predict the evolution of a conserved scalar variable in a 
statistically homogeneous field, in a manner that is in accord with both DNS and 
experimental measurements. Despite this positive feature, some questions 
pertaining to the generality of this density for more realistic flow configurations 
remained unanswered. 

In this work, we intend to address some of these remaining questions. 

Specifically, we shall make a comparative assessment of the behavior of this 
model, and provide an estimate of the degree of its validity in the context of more 
general conditions than previously considered. For this purpose, we have selected 
a configuration of particular interest in the design of diffusion controlled reactors; 
namely, a spatially developing parallel mixing layer under the influence of 
harmonic perturbations. The flow field which develops in this setting is dominated 
by large scale coherent structures and exhibits non-homogeneous features. This 
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inhomogeneity plays a dominant role in the mechanism of mixing in turbulent 
flows and is very suitable for promoting mixing and chemical reactions. Because 
of this property, this configuration has been widely utilized in numerous 
investigations on turbulent reacting flow phenomena (see Givi and Riley (1992) 
for a recent review). 

We shall follow an analogous approach to that of our previous work (Madnia et 
al., 1991a) in assessing the model’s behavior. However, there are fundamental 
differences between the works in formulating the problem and in implementing 
the numerical procedure involved in the DNS. The configuration considered here 
is non-homogeneous. Therefore, both the procedure of statistical sampling and 
the numerical algorithm used in the DNS are different from the homogeneous 
simulations of Madnia et al. Also, the mathematical operations leading to the 
estimate of the reactant conversion rate are substantially more elaborate. 
However, analytical solutions can be attained with the aid of the mathematical 
properties of the Beta Function and the Incomplete Beta Function. With these, 
the generated results can be directly compared with those of DNS. 


DESCRIPTION 


The problem under consideration is that of a spatially evolving mixing layer 
containing initially segregated reacting species A and B in the two free streams. 
Species A is introduced into the upper, high-speed stream and species B enters on 
the lower, low-speed side. The chemical reaction between the two species is 
assumed to be single step, irreversible and infinitely fast. The magnitude of the 
Mach number is assumed negligible. Therefore, the flow is considered incompres- 
sible. The two species are assumed thermodynamically similar with identical 
diffusion coefficients. Within the framework of these approximations, the flow is 
mathematically described by the conservation equations of mass, momentum, and 
a species conservation equation for the trace of a conserved Shvab-Zeldovich 
variable J?, which characterizes the compositional structure of the flow. These 
equations are parameterized by the Reynolds number, the Peclet number, and 
the velocity ratio across the layer. 

The computational package employed in the DNS is based on a hybrid 
pseudospectral-spectral element algorithm developed by Givi and Jou (1988) and 
McMurtry and Givi (1991). The pseudospectral routine, which is based on 
Fourier collocation, is used in the cross stream direction together with free slip 
boundary conditions. The spectral-element discretization is employed in the 
streamwise direction. This involves the decomposition of the domain in this 
direction into an ensemble of macro finite elements. Within each of these 
elements, the dependent flow variables are approximated spectrally by means of 
Tchebysheff polynomials of the first kind. The hybrid procedure employed in this 
way is very attractive in that it combines the accuracy of spectral discretization 
with the versatility of finite element methods. Therefore, it is convenient for 
accurate simulations of the complex spatially developing flow under consideration 



MODELING OF REACTING TURBULENT SHEAR FLOWS 199 

here. For a detailed description of this hybrid method for DNS of turbulent shear 
flows, the reader is referred to McMurtry and Givi (1991). 

The flow field is initialized at the inflow by a hyperbolic tangent mean velocity 
distribution, upon which low amplitude perturbations are superimposed in order 
to promote the formation of coherent vortices. The initialization of the 
Shvab-Zeldovich variable at the inflow is such that it takes on the values of 
zero and one in the streams containing B and A, respectively. The disturbances 
on the mean flow correspond to the most unstable mode of the hyperbolic tangent 
profile (Michalke, 1965) and its first two subharmonics. The amplitude of the 
fluctuating velocity is set equal to 6% of the mean velocity. This amplitude 
corresponds to that measured in typical shear layer experiments. The magnitudes 
of the phase shifts between the modes of the instability waves are randomly 
selected from a random seed with a top hat PDF of zero mean and a specified 
variance. The implementation of these phase shifts is the only mechanism to 
introduce randomness into an otherwise deterministic simulation. This is to 
partially mimic the random “commotions” which are present in the universe, 
into an isolated two-dimensional simulation. At the outflow, a weak condition of 
zero second derivative is applied for all the dependent variables. This boundary 
condition cannot be justified in a rigorous mathematical sense and can be 
specified only in an ad hoc manner. This was done in order to facilitate 
implementation of the Dirichlet boundary conditions in the finite element 
procedure employed here. Moreover, the results of extensive numerical tests 
showed that the effects of the approximate boundary conditions are confined 
within the last computational element. This is consistent with that found 
previously by Korczak and Hu (1987). Therefore, the solution in the last 
computational domain can be ignored. 

With the solution of the unsteady transport equations, the DNS data are 
extremely useful in visualizing the instantaneous behavior of the flow. The results 
of these direct simulations can also provide an assessment of the statistical 
behavior of the flow. This is due to the availability of flow information at all the 
computational grid points and at all times. Therefore, statistical sampling can be 
done quite easily by ensemble averaging of the instantaneous data gathered over 
many realizations. 

In our simulations the instantaneous data of the Shvab-Zeldovich variable, 
obtained from the DNS, are used for statistical analysis. A total of 3200 
realizations were employed in the sampling. With this ensemble, the magnitudes 
of the mean and the variance of the Shvab-Zeldovich variable are calculated at 
all the grid points. These statistical quantities are of primary interest for 
turbulence modeling. Based on the knowledge of these first two moments, a Beta 
density is assumed for the PDF of the Shvab-Zeldovich variable. The ensemble 
mean values of the product concentration and the unmixedness are calculated 
subsequently via this density. These modeled quantities are then compared with 
those obtained directly from the DNS data. At the majority of grid points the 
Beta density is asymmetric. Nevertheless, analytical solutions for the statistics of 
the reacting field are possible and have recently been obtained by Madnia et at. 
(1991b). The derivations of the equations are rather involved. Here, we only 
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present the final results for the normalized mean product concentration (C p ) and 
that of the unmixedness W 2 : 


(Cp) = (C p )(x,y) = l-e l -6 2 (1) 

'V 2 = 'V 2 (x,y) = -6 l d 2 (2) 


where: 


1 r(tv + p) a~P 
1 2 a+fi ~\a + 0) r(a)T(P) + or + 0 (l ~ I ^° C ' 


B,= 


f (a + P) p - a 


2 a+p ~ l (a + P) r(a)r(P) + a + p ^ 


(3) 

(4) 


Here, T is the Gamma function, and a and P are related to the first two moments 
of the random field (Frankel et al, 1991). /*.(<*, 0) is the Incomplete Beta 
Function, and / s , denotes the stoichiometric value of the Shvab-Zeldovich 
variable, /. For unity normalized concentrations at the free streams, / sl = 0.5. 
The results predicted by these equations shall be compared with those generated 
by DNS for a quantitative assessment of the closure. 


PRESENT ATATION OF RESULTS 

Computations were performed in a domain of size 13<5 x 34<5, where 6 denotes 
the vorticity thickness at the inlet of the flow. There are 65 Fourier modes in the 
cross-stream direction and 42 finite elements are used for the streamwise 
discretization. A fifth order Tchebysheff polynomial is employed to approximate 
the variables within each element. This discretization is equivalent to, at least, a 
fifth order accurate finite difference technique even if the spectral convergence is 
not considered. This results in a total number of 211 points in the streamwise 
direction. A second-order Adams-Bashforth finite difference scheme was utilized 
for temporal discretization. With the computationally affordable 211x65 grid 
resolution available, accurate simulations with Re = Pe = 200 (where these 
normalized quantities are defined based on the velocity difference across the 
layer, and the initial vorticity thickness) are possible (Givi, 1989). These values 
are somewhat smaller than those of a fully developed laboratory turbulent shear 
layer, but are sufficiently high to promote the growth of instability waves. 

In order to visualize the flow evolution, a time series of DNS generated 
snapshots of the instantaneously normalized product concentration are shown in 
Figure 1. This figure depicts the way in which the small amplitude perturbations 
manifest themselves in the formation of large scale coherent vortices downstream 
of the splitter plate. The forcing associated with the most unstable mode alone 
would cause the initial roll-up of these vortices. The perturbations corresponding 
to the first subharmonic mode result in a second roll-up in the form of merging 
neighboring vortices, and the presence of the second subharmonic generates a 
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FIGURE 1 Time sequence of contour plots of instantaneous product concentration. 

(See color Plate I.) 

third roll-up (second merging) at a region near the outflow. The roll-up and the 
subsequent pairings of these vortical structures engulf the unreacted species from 
the two streams into the mixing zone. The dynamics of the vortical evolution in 
this process reflect the two main mechanisms of mixing enhancement. Large 
concentrations of vorticity, that are approximately of one sign, bring the fluid 
from the two free streams into the large scale structures. Within these structures, 
the fluid elements are subject to further straining, and diffusion between the two 
fluids results in final mixing at the molecular scale. With these processes, along 
with the assumption of infinitely fast chemistry, the rate of product formation is 
naturally enhanced by the dynamics of vortical evolution. 

The ensemble average of the product concentration, formed as a result of the 
mechanisms just described, is presented in Figure 2(a). The statistical averaging 
process was performed over a time span approximately equal to the residence 
time of the flow within the domain of computational consideration. This analysis 
was performed after the effects of the initial transients were washed away at the 
outflow. A comparison between Figures 1 and 2(a) reveals the essential 
deficiencies of the ensemble averaged results. Note that, all the detailed 
intricacies of the dynamics of flow evolution are masked by the averaging. While 
these ensemble-averaged results are of main interest for practical applications, 
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FIGURE 2 Contour plots of mean product concentration generated by (a) DNS, (b) the model 

(See color Plate II.) 

they are not capable of describing many of the interesting features of turbulence 
transport. The results depicted in Figure 2(a) are the best that one can expect to 
generate from a turbulence closure in a statistical sense. 

To demonstrate the capability of the model based on the Beta density, in 
Figure 2(b), the average product concentration obtained by means of Eqs. 
(1)_(4) is presented. A comparison of this figure with Figure 2(a) reveals a 
remarkable similarity between the two results. For a more quantitative com- 
parison, the cross stream variation of this product concentration and those of the 
negative of the unmixedness are presented at two streamwise locations in Figures 
3-4. Also the streamwise variation of the total product, T p , is presented in Figure 
5. In all these figures, the comparisons between the model predictions and the 
DNS results are noteworthy. 

The agreements demonstrated by the comparisons noted above provide a 
reasonable justification for recommending Eqs. (1)— (4) as working relations in 
routine engineering predictions. However, these relations are advocated only for 
statistical predictions of low order moments. The physics of the mixing phenome- 
non in a spatially developing shear flow is far too complex to be completely 
described by the Beta density. The highly intermittent nature of the mixing 
process exhibits features that cannot be fully captured by this density (Givi, 
1989). To demonstrate this point quantitatively, in Figures 6-7 the cross stream 
variations of the third and fourth centralized moments of the Shvab Zeldovich 
variable calculated from the Beta density are compared with those generated by 
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Collocation points 

FIGURE 5 Streamwise variation of total product (7^) from the DNS and the model. 

DNS. These figures indicate that while the trend is somewhat similar in the two 
cases, the extent of agreement is not as good as those of lower level statistical 
quantities (Figures 2-5). This is primarily due to the mechanism of large scale 
entrainment in the shear flow, in that there are always unmixed fluids present 
within the core of the vortical structures. This results in a slight trimodal behavior 
in the probability distributions generated by DNS (Givi and Jou, 1988; Lowery 
and Reynolds, 1990) and observed experimentally (Koochesfahani, 1984, Dimota- 
kis, 1989). This trimodal behavior cannot be captured by the Beta density which 
is at best suitable for producing bimodal distributions. 

Recall that the specification of the Beta density is based on the knowledge of 
the first two moments of the conserved scalar. Here, these moments were 
supplied via DNS, but in an actual engineering application they must be provided 
by other means ( i.e . an appropriate turbulence closure, or experimental data). 
However, specification of these parameters is substantially simpler than those of 
the reacting scalars, even with the assumption of an infinitely fast reaction. This is 
simply due to the fact that such data can be obtained in simpler, non-reacting flow 
experiments without any physical complications associated with the chemistry. 
With such data available, then Eqs. (1)— (4) are deemed very suitable, at least in 
the absence of better alternatives, for a reasonably accurate and inexpensive 
prediction of the limiting bounds of the reactant conversion rate. These equations 
are therefore recommended for this purpose in either homogeneous or non- 
homogeneous flows. 
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FIGURE 7 Cross stream variations of fourth moment (ji 4 ) of the Shvab-Zeldovich variable from 
the DNS and the model; (a) x /d = 11, (b) x/d = 22. 
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NOMENCLATURE 


A,B Reactants. 

C p Product Concentration. 

$ Shvab-Zeldovich variable. 

Pe The Peclet number. 

Re The Reynolds number. 

T p (x) Total product = jlZdy $odx(C p )(x, y). 

x,y The physical coordinates. 

<5 Vorticity thickness at the inflow. 

W 2 The unmixedness. 

H 3 Third order moments of the Shvab-Zeldovich variable. 

/u 4 Fourth order moments of the Shvab-Zeldovich variable. 

() Ensemble average. 
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The Compositional Structure and the Effects of Exothermicity 
in a Nonpremixed Planar Jet Flame 

C. J. STEINBERGER,* T. J. VIDONI, and P. GIVI 

Department of Mechanical anil Aerospace Engineering. State University of New York at Buffalo , 

Buffalo, NY 14260-4400 

Results are presented of direct numerical simulation (DNS) of a randomly perturbed compressible, spatially 
developing two-dimensional planar jet under the influence of a finite rate chemical reaction of the type 
p + o - Product. The objectives of the simulations arc to assess the compositional structure of the flame 
and to determine the influence of reaction exothermicity by means of statistical sampling ol the 
DNS generated data. It is shown that even with this idealized kinetics model the simulated results exhibit 
features in accord with experimental data. These results indicate that the Damkohler number is an important 
parameter in determining the statistical composition of the reacting field and that the results are not very 
sensitive to the mechanism by which this parameter is varied. It is demonstrated that as the intensity of 
mixing is increased and the effect of finite rate chemistry is more pronounced, the magnitudes of the 
ensemble mean and variance of the product mass fraction decrease and those of the reactants mass fraction 
increase. Also, at higher mixing rates the joint probability density functions of the reactants mass fractions 
shift towards higher values within the composition domain, indicating a lower reactedness. These trends are 
consistent with those observed experimentally and arc useful in portraying the statistical structure of 
nonequilibrium diffusion flames. The DNS-generated data are also utilized to examine the applicability of 
the "laminar diffusion flamelet model” in predicting the rate of the reactant conversion with finite rate 
chemistry. This examination indicates that the performance of the model is improved as the value of the 
Damkohler number is increased. Finally, the simulated results suggest that in the setting of a "turbulent" 
flame, the effect of the heat liberated by the chemical reaction is to increase the rate of reactant conversion. 
This finding is different from those of earlier DNS results and laboratory investigations that indicate a 
suppressed chemical reaction with increasing heat release. 


INTRODUCTION 

Recent advances in flow visualization and diag- 
nostic techniques have made significant con- 
tributions to the understanding of turbulent 
combustion phenomena [1, 2]. With imple- 
mentation of these advanced techniques it is 
now possible to make a detailed examination 
of the structure of turbulent flames and to 
establish a better understanding of the intri- 
cate physics of such flames. Among recent 
applications of these diagnostic techniques, 
Masri et al. [3, 4] have made noteworthy 
progress in describing the structure of non- 
premixed turbulent jet flames under the influ- 
ence of finite rate chemical reactions (see Refs. 
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5 and 6 for comprehensive reviews). In particu- 
lar, they have reported results provided by 
statistical sampling of data obtained by 
Raman-Rayleigh scattering measurements of 
the species mass fractions and the temperature 
in a turbulent methane jet flame. These mea- 
surements are utilized primarily for construct- 
ing single-point statistics of the reacting field 
including the ensemble mean, the variance and 
probability density functions (pdfs) of the rele- 
vant thermochemical variables. The response 
of the flow under different mixing conditions is 
characterized by a comparison of the results 
obtained at various mixing rates. The sampling 
of the data in this way allows a systematic 
assessment of the compositional structure of 
the nonequilibrium flame. Such data are valu- 
able in describing the phenomena of mixing 
and chemical reactions in nonpremixed turbu- 
lent combustion and in addressing some of the 
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existing problems in the modeling of such 
flames oy means of statistical methods. 

In order to provide a computational comple- 
ment to laboratory investigations, there has 
been significant progress in the development 
and implementation of direct numerical simu- 
lation (DNS) of nonpremixed reacting flows 
[7-14] (see Refs. 15-19 for recent reviews). 
The extent of progress made by these contribu- 
tions has been very encouraging, thus justify- 
ing further utilization of direct methods in 
the analysis of turbulent flames. In this work 
we intend to continue our investigation of 
nonpremixed reacting flows via DNS and to 
simulate such flows in a context relevant to 
laboratory investigations. Obviously, it is still 
impossible to perform DNS of a realistic 
“physical” flame [15]. Here, our intention is to 
make use of present computational capabilities 
to further address some of the fundamental 
and pertinent issues in regard to the dynamics 
of these flows. 

In the next section the method of investiga- 
tion is presented together with the underlying 
assumptions imposed in making the problem 
tractable for computational treatment. The 
ramifications associated with some of the 
assumptions are discussed, where appropriate, 
in the results. The focus is on assessing the 
compositional structure of the flame and on 
revealing the influences of reaction exother- 
micity. It is shown that some of the trends 
observed in laboratory experiments can be 
captured by DNS, and speculation is made 
regarding some of the other features not yet 
ascertained in laboratory investigations. 

METHOD OF INVESTIGATION 

The flow under investigation consists of a spa- 
tially developing, two-dimensional, planar jet 
flame. The schematic diagram of the jet is 
shown in Fig. la along with the specifications 
of physical dimensions. The fuel, F, discharges 
from the inner higher speed jet of width D 
with a velocity of U F into a co-flowing lower 
velocity (U 0 ) oxidizer, O, stream. The hydrody- 
namic field is characterized by the velocity 
ratio r = U F /U 0 , the Reynolds number (Re) 
based on the velocity difference, A(7 = U F — 
U Q , across the layer (Re = pAt/D//x), and the 


convective Mach number in each stream. Tin 
planar configuration considered here is diffei 
ent from that of a circular jet flame in Refs, 
and 4, but allows the elucidation of some spe 
cific physical features that are believed to lx 
invariant of the jet geometry [18]. The flow 
evolves spatially in the streamwise direction 
(.v), and impermeable free slip boundary con 
ditions are imposed in the cross-stream direc- 
tion (y). The DNS computational procedure is 
the same as that employed in Refs. 14 and 20. 
The numerical discretization procedure is 
based on a two-four (second-order accurate in 
time and fourth-order truncation in space) 
compact parameter finite difference scheme. 
This scheme has been developed and made 
operational by Carpenter [21], and is utilized in 
the context of the SPARK computer code [22]. 
A generalized single-step chemistry model of 
the form F + O -» Product is assumed. This 
model is simple enough to allow an efficient 
treatment via DNS and yet it is capable of 
capturing some of the important nonequilib- 
rium effects associated with combustion. 
All the species involved in this reaction are 
assumed to be thermodynamically identical and 
the values of the physical transport parameters 
are assumed to be invariant with temperature. 

The flowfield is initialized with a top hat 
mean streamwise velocity profile with a low 
amplitude random forcing at the inner jet 
entrance ( x = 0). This forcing is imposed across 
the layer and is in terms of sinusoidal distur- 
bances of specified frequency and phase shifts. 
The magnitudes of the frequency' and the phase 
shifts between the different forcing modes are 
selected randomly from a seed with a uniform 
probability distribution and with specified val- 
ues of the mean and variance. The random 
specification of the frequency and the phase 
shifts is intended to mimic the influence of 
turbulence in an otherwise deterministic simu- 
lation. The imposition of external perturba- 
tions expedites the formation of large scale 
structures; no attempt is made to investigate 
stability characteristics of the jet. The ampli- 
tudes of the perturbation are specified such 
that the maximum intensity of fluctuation at 
the inflow is no more than 1.5%. The inlet 
temperature ( T ) and pressure ( P ) are the same 
in both streams, and the reactants are com- 
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ORIGINAL FA Of is 
OF POOR QUALITY 


pletely segregated at the inlet. To initiate reac- 
tion, the temperature at the interface between 
the two reactants is increased at the inlet. An 
approximate Gaussian profile is adopted tor 
the temperature distribution, and an error 


function profile is assumed for the mass trac- 
tions at the interfaces of the two reactants. 
With the assumption of a second-order chem- 
istry model the rate of reactant conversion ( a> ) 
is expressed by to = K, C, C (r Here, C and K, 
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Fig. 1. (a) Schematic diagram of the planar jet configuration, (b) The grid. 
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denote the concentration and the reaction rate 
constant, respectively. Both constant rate 
kinetics and an Arrhenius chemistry model are 
considered. In the former K f is constant, and 
in the latter K t = A f exp( — E/y?T). In this 
expression denotes the preexponential 

factor, E is the activation energy, and 
represents the universal gas constant. The 
nondimensional parameters influencing the 
rate of reactant conversion in this setting 
are the Damkohler number [Da = K f C F / 
(1U/D)], the Zeldovich number in the 
Arrhenius model (Ze = E/^T t ), the heat 
release parameter (Q = -AH°/C y T F ), the 
Prandtl (Pr) number, and the Schmidt (Sc) 
number. In these expressions, the subscript i 
denotes the inlet, C v is the specific heat at 
constant volume, and -A// n indicates the 
enthalpy of combustion. 

Several simulations are conducted. Table 1 
provides a listing of the values of the physi- 
cal parameters considered in each simula- 
tion. In these simulations, the influences of 
the Reynolds number, the velocity ratio, the 
Prandtl number, and the Schmidt number are 
not investigated. The values of the molecular 
Prandtl and the Schmidt numbers are set to 
unity (Pr = Sc = 1) to make the analysis asso- 
ciated with laminar diffusion flamelet model- 
ing easier. The magnitude of the Zeldovich 
number is held constant in all the Arrhenius 
simulations but the values of the Damkohler 
numbers and the heat release parameter are 
altered to assess their influence on the dynam- 
ics of combustion. With the assumption noted 
above, the speed of sound ( a ) is uniform at the 
inflow. Therefore, the magnitude of the con- 

TABLE 1 


Physical Parameters for Simulations 


Flame 

K 

Da 

Ze 

Q 

A 

0.1 

0.1 

20 

0.5 

B 

0.1 

0.05 

20 

0.5 

C 

0.2 

0.05 

20 

0.5 

D 

0.1 

0.05 

0 

0.0 

E 

0.1 

0.5 

0 

0.0 

F 

0.1 

10.0 

0 

0.0 

G 

0.1 

100.0 

0 

0.0 

H 

0.1 

0.1 

20 

0.0 

I 

0.1 

0.1 

20 

1.5 

J 

0.1 

0.1 

20 

6.0 


vective Mach number is the same in both 
streams at the inlet and is given by ,1/ = (E 
U I )/ a where the convective velocity, is 
expressed as U c = (U F + U G )/ 2. The full 
compressible form of the transport equations 
[23] are considered in the computational for- 
mulation. However, the magnitudes of the con- 
vective Mach numbers considered are not very 
large. Therefore, many of the physical issues 
associated with high-speed combustion [14, IS, 
22] are not addressed. Also, the body forces 
are assumed negligible; therefore, the inter- 
esting physical problems associated with the 
buoyancy [24, 25] are not considered. 

The compositional structure of the flame is 
assessed by means of examining the response 
of the flow to the Damkohler number in a 
manner similar to that in Ref. 3. This assess- 
ment is quantified by considering the statistical 
properties of the reacting field including the 
ensemble mean, the variance, and the pdfs 
of the mass fractions. The analysis has been 
made in a heat releasing reacting jet with an 
Arrhenius chemistry model. Flames A, B, and 
C (Table 1) are considered for this purpose. 
These flames are identified by the relative 
magnitude of the mixing intensity quantified 
by the value of the Damkohler number. In 
Flames B and C, the magnitudes of the 
Damkohler number are the same and are set 
equal to half of that in Flame A. However, the 
mechanism by which this is enacted is differ- 
ent. In Flame B, this is implemented by 
decreasing the magnitude of the chemical 
reaction frequency, and in Flame C by increas- 
ing the magnitude of the characteristic velocity 
of the jet. Both of these procedures result 
in the same reduction in the value of the 
Damkohler number. The magnitude of 
the convective Mach number is higher in Flame 
C than that in Flame A, but the value of this 
Mach number is not so large as to have a 
substantial influence on the spatial growth of 
the jet [14]. Also, the imposition of external 
forcing masks the effects of other possible mix- 
ing inhibition mechanisms (such as exothermic- 
ity, density stratification, etc.) [12, 14, 26-29], 

To examine the performance of the diffusion 
flamelet model and to assess the effects of heat 
release. Flames D through J are considered 
(Table 1). In Flame D through G, a constant 
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rate non-heat-releasing chemical reaction 
is simulated. In Flames H through J, an 
Arrhenius reaction is considered with varying 
magnitude of heat release. 

RESULTS 

The simulations are performed on the mesh 
shown in Fig. lb. The grid generation proce- 
dure involves a nonuniform mesh with a heavy 
concentration of grid points near the regions 
of maximum instantaneous shear [14, 20, 22]. 
With the computational resources available to 
us, routine numerical experiments with a reso- 
lution of 245 x 165 grid points are feasible. In 
some of the runs the independency of the 
results from the grid spacing was assessed by 
performing simulations with 405 x 255 grid 
points. This assessment was made for Flames 
A through C and the results indicate the con- 
vergence of the statistical results. Therefore, 
the simulated results presented here are based 
on a resolution consisting of 405 X 255 grids 
for Flames A through C, and of 245 X 165 
grids for Flames D through J. These simula- 
tions are conducted with moderate values 
of the Reynolds, Peclet, Damkohler, and 
Zeldovich numbers within a domain L x = 
13 \D > jc > 0, L y = 6f D > y > 0. With this 
resolution it was possible to perform simula- 
tions with r = 2 and Re = 1.0 X 10 4 free of 
spurious oscillations. The smallest grid spacing 
in the cross stream direction is near the center 
of the layer and is Ay = 6 X 10 -5 m. This 
spacing is smaller than the diffusion length 
scale. However, at regions far from the center 
of the layer the spacing is relatively coarse. At 
these regions obviously the scales are not fully 
resolved, but the magnitudes of the gradients 
here are not that large to cause numerical 
oscillations. This is shown in the next section. 
The computational time increment varies at 
each time step and its value is the smallest of 
the characteristic time scales associated with 
hydrodynamics and chemistry. As shown in 
Table 1, the magnitudes of the Damkohler 
numbers are small enough to avoid the need 
for prohibitively small time increments. With 
the resolution adopted, a typical simulation 
requires about 35 h of CPU time on a Cray-2 
supercomputer. 


To visualize the global structures, contour 
plots of the instantaneous fuel mass fractions 
arc shown in Fig. 2a. This figure shows how the 
growth of low-amplitude random perturbations 
manifest themselves in the formation of large 
scale coherent vortices downstream of the jet 
entrance. Due to the random nature of the 
forcing mechanism, the evolution of the large- 
scale structures is complex and, similar to 
laboratory shear flows, consists of random 
formation of vorticity rollups and subsequent 
pairing and coalescence of neighboring vor- 
tices. Also, due to the response of the layer to 
sinuous instabilities the layer is asymmetric 
with respect to the jet centerline. The forma- 
tion of large-scale structures results in an 
enhancement of mixing. Large concentrations 
of vorticity on both sides of the jet centerline 
bring the fluid from the two streams into the 
core of the large scale structures. Within these 
structures the fluid elements are subject to 
further straining and the molecular diffusion 
between the two fluids results in final mixing at 
the molecular scale. The region of large gradi- 
ents is confined to the inner jet zone where 
there is heavy concentration of grid points and 
there are no numerical oscillations in resolving 
these gradients. The enhanced mixing mecha- 
nism caused by large-scale vortices results in a 
general depletion of fuel. Near the entrance of 
the jet where the reactants are first brought 
into contact, the extent of chemical reaction is 
fairly uniform at the reactants’ interface. Fur- 
ther downstream, as the magnitude of the in- 
stantaneous strain increases, the reaction rate 
approaches zero at the braids and the flame 
locally quenches at those locations. This mech- 
anism of flame extinction is consistent with 
previous DNS results [9, 10], and is also cor- 
roborated by experimental observations [30, 
31], This behavior was observed in all Arrhe- 
nius simulations with moderate Damkohler 
numbers indicating the departure from chemi- 
cal equilibrium. 

Compositional Structure 

To examine the effects of finite rate chemistry 
on the compositional structure, the results 
extracted by DNS are statistically analyzed. 
With the assumption of statistically stationary 



Fig. 2. (a) Plot of product mass fraction contours, (a) Instantaneous data, (b) Ensemble-averaged data. 
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(low, the instantaneous data can be used to 
construct ensemble-averages. Flere, these aver- 
ages arc obtained by statistical analysis of 9000 
samples gathered after the effects of initial 
transients are washed away at the outflow. The 
data sampling was performed within a period 
equivalent to 3 times the residence time of the 
flow. With this sampling, the maximum errors 
in evaluating the mean and the variance of the 
scalar quantities are estimated to be 4.3% and 
8.7%, respectively. Typical results of such aver- 
aged quantities are presented in Fig. 2b in the 
form of a contour plot of ensemble-averaged 
fuel mass fraction. 

To exhibit the effects of finite rate chemistry 
in a quantitative manner, comparisons are 
made between the statistical data in Flames A, 
B, and C. In Fig. 3 the cross stream variation 
of the ensemble-mean values (denoted by ( )) 
of the product mass fraction is shown at a 
streamwise location (x = 10D). Figure 3 shows 
that the magnitude of mean product mass frac- 
tion in Flame A is higher than those of Flames 


B and C. Phis suggests that as the llame 
approaches extinction, the mass fractions ot 
the reactants in the reaction zone increase and 
the rate of product formation decreases. TIk 
same behavior is observed at all streamwise 
locations, as indicated by the spatial variation 
of the integrated product concentration T ( y . 
shown in Fig. 4. Here T Cp = f(C,,)ciy. A simi- 
lar behavior is observed in Fig. 5, which shows 
that the amplitude of the variance of the prod- 
uct mass fraction is higher in Flame A than 
those in Flames B and C. The trend observed 
in Figs. 3-5 is consistent with that in the 
laboratory measurements [3] and reveal the 
effects of the mixing intensity on compositional 
structure of the flame. These results also sug- 
gest that this behavior is somewhat insensitive 
of the procedure by which the Damkohler 
number is varied. In the experiments, the 
Damkohler number was altered only by means 
of increasing the hydrodynamic intensity (i.e., 
increasing the inner jet velocity). This was 
because a single fuel was used and, therefore. 



224 


C. J. STEINBERGER ET \I 


the chemical frequency was constant in all the 
experiments. Our results indicate that this 
dependence is observed regardless of the pro- 
cedure by which the change is made and Flames 
B and C portray approximately the same statis- 
tical conduct. The behavior observed in Figs. 
3-5 demonstrates the dominant influence of 
the Damkbhler number in quantitative analysis 
of nonequilibrium flames and also in the com- 
parative assessment of the statistical results 
generated in our DNS with those obtained in 
laboratory experiments (See Ref. 32 for addi- 
tional discussion concerning this issues). 

The influence of finite rate chemistry on the 
compositional structure of the flame is further 
considered by examining the marginal and the 
joint pdfs of the reacting scalar variables. In 
Fig. 6, results are presented of typical single- 
point joint pdfs of the fuel and the oxidizer 
mass fractions, i.e., £P(Y F ,Y 0 ) for Flames A 
and B. At a lower mixing rate. Flame A (Fig. 
6a), there is a greater amount of product 


formed and the pdfs are higher at the mixed 
mass fraction values at the interior ol 
the allowable region (0 < Y h , Y 0 < 1). At the 
higher mixing rate, Flame B (Fig. 6b), the 
mixture is composed of the two reactants with 
the joint pdfs concentrated closer to the line of 
Y h + Y a = 1. These results are consistent with 
those reported in Ref. 4 and indicate that at 
the higher levels of mixing intensity the results 
are closer to those obtained by cold mixing. 
The corresponding single-point marginal pdfs 
of the product mass fraction ^(Vj,) are shown 
in Fig. 7 for Flames A, B, and C. In all of these 
flames the pdfs are concentrated near zero at 
free streams and become broader near the 
central region of the layer. As the mixing rate 
is decreased the pdfs are shifted toward higher 
values of mass fraction. Again, the pdfs for 
Flames B and C are very similar. The trend in 
Fig. 7 is consistent with the experimental data 
[4] and indicates lower reactedness at higher 
muting rates. However, due to strong intermit- 
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Fig. 4. Streamwise variation of the total product concentration. 
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tency of the flow in the two-dimensional simu- 
lations the pdfs do not adopt a Gaussian distri- 
bution, and are heavily populated near free 
stream values, that is, zero. Also, the tails of 
the pdfs become thinner as the mixing is inten- 
sified. This suggests a reduction in the variance 
of the product mass fraction as was reported in 
the experiments of Ref. 3 and also evident in 
Fig. 5. 

Laminar Diffusion Flamelet Model 

The data bank generated by DNS provides an 
effective means of assessing the performance 
of turbulence closures, at least in the setting of 
simple flows. In fact, within the past decade 
DNS-generated data have proven valuable in 
such assessments [15, 33]. One of the most 
recent closures for the analysis of nonequilib- 
rium flames is the “laminar diffusion flamelet 
model” proposed in Refs. 31, 34, and 35. 
According to this model, a turbulent flame is 
viewed as an ensemble of thin laminar diffu- 


sion flamelets. With this view, effectively, the 
influences of hydrodynamics and chemical 
reactions are decoupled and the compositional 
structure of the flame is assumed to be the 
same as that in a laminar flame configuration 
with the same chemical characteristics. In the 
limit of an infinitely fast chemistry, or chemical 
equilibrium, this is obvious since the composi- 
tional structure depends on the magnitude of 
the mixture fraction. That is, the instantaneous 
values of a reacting scalar ip is uniquely deter- 
mined from the local values of the mixture 
fraction f, i.e., ip = i K£(x,/» [23], In finite 
chemistry, the laminar flamelet model implies 
that nonequilibrium effects manifest them- 
selves through a dependence on a second 
parameter, namely, the dissipation rate of the 
mixture fraction, x ■ That is, ip ~ X^- 

In order to examine the performance of the 
model in describing the compositional struc- 
ture of the jet flame, simulations are per- 
formed of Flames D through G (Table 1). For 
ease of analysis in these simulations the reac- 
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tion rate (Kj ) is assumed constant but is dif- 
ferent in each flame, and the magnitude of 
heat release is set equal to zero. For the pur- 
pose of comparison with the flamelet model, a 
simple laminar system consisting of a steady 
one-dimensional opposed jet laminar flame 
configuration [36-38] is considered. In this set- 
ting, with the assumption of unity Lewis num- 
ber the simulations based on the flamelet 
model involves the solution of the diffusion- 
reaction equation [38]: 

d^_ - 2 ^ 

d£ 2 PX ' ^ 

Here, denotes the chemical source term 

corresponding to scalar ip, and the variable x 
is related to £ through 

* = * s . exp[ — 2(erT , (2£- 1))~] . (2) 


where ;y st denotes the stoichiometric value of 
X- In the laminar flame, this stoichiometric 
value is the same as the “stagnation value” of 
X and is determined by the hydrodynamic in- 
tensity of the jet. However, for the unsteady jet 
flame considered in DNS, a single value cannot 
be adopted. Instead, only some representative 
values can be considered. With the instanta- 
neous values of (£, *) extracted from the DNS 
results, the flamelet library is constructed by 
means of Eqs. 1-2. With this construction the 
results are described in terms of </r(£, * st ). In 
the laminar flame, these results are obviously 
single valued functions of the local Damkohler 
number, K^/ x it , as indicated by Eqs. 1—2. In 
the unsteady jet flame, the results are extracted 
directly by DNS and are parameterized in terms 
of the instantaneous values of K,/x sX . The 
comparison between the model and DNS 
results for one of these flames is shown in Fig. 
8, which presents the “scatter plot” of product 
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mass fractions versus the mixture fraction. For 
the laminar flame configuration, each of the 
curves with a constant value of K f /\ st repre- 
sents a flamelet. For DNS, the compositional 
structure is constructed by parameterization of 
the data in terms of the instantaneous values 
of x ■ If the laminar diffusion flamelet model is 
to be considered valid, then all the scatter data 
in the range C, < K f /x %t < C 2 should fall be- 
tween the laminar flamelet solution corre- 
sponding to K f /x u = C, and K f / Xtt = C 2 . 
Figure 8 indicates that at low values of the 
local Damkohler number the model underpre- 
dicts the DNS data. The agreement becomes 
somewhat better at high values of this number. 
However, in all ranges of K f /x st , there is a 
substantial portion of DNS data that fall out- 
side the boundaries suggested by the flamelet 
model. The results are bounded in the upper 
limit as the value of the local Damkohler num- 
ber approaches infinity. In this flame sheet 
limit both results obviously coincide and are 
single-valued functions of £. A more clear 


comparison is made in Fig. 9, which shows the 
scatter of product mass fractions versus the 
local Damkohler number, within the window of 
0.45 < £ < 0.55, that is, £ = g s = 0.5. Figure 9 
clearly shows that as the Damkohler number is 
increased, the data points shift towards the 
right and fall between the window of the two 
solutions predicted by the flamelet model. 
Based on these results, it is concluded that the 
laminar flamelet model predicts the composi- 
tional structure of the flame accurately at high 
Damkohler numbers, that is, near the limit of 
an infinitely thin flame. The use of the model 
in flames with strong departure from this limit 
cannot be guaranteed for general applications. 

Exothermicity Effects 

Despite the attractive features of our simula- 
tions in predicting the trends observed experi- 
mentally, it must be emphasized that in the 
context considered the factors influencing the 
stability characteristics of the jet are masked 
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Fig. 7. Marginal pdfs of product mass fraction, Y P ) at x = 10 D, y = 3.8 D. 


by addition of external forcing [28], Note that 
these simulations are of a transitional flow, 
whereas the experimental data are obtained in 
a fully turbulent jet flame. The imposition of 
explicit forcing was therefore necessary to 
mimic a turbulent environment. This forcing, 
however, overcomes the mixing suppression 
caused by factors such as reaction exothermic- 
ity, compressibility, etc. [14, 26-29]. In order to 
make this point clear, simulations are per- 
formed of several heat releasing flames (Flames 
H through J in Table 1). In all these flames, 
the magnitudes of the velocities, the 
Damkohler number and the Zeldovich number 
are identical and are set equal to those in 
Flame A. However, the extent of exothermicity 
is varied. The total integrated product mass 
fraction T V/> (x) generated by these flames is 
compared in Fig. 10, which reveals a higher 
rate of product formation as the magnitude of 
the heat release is increased. This observation 
is not consistent with any of the previous DNS 


[8, 14] or experimental [39] results. The reason 
for this discrepancy is due to the effects of the 
Arrhenius kinetics model in the present simu- 
lations and cannot be observed in constant 
rate kinetics simulations and/or experiments 
[8, 14, 39], The increased level of heat release 
also delays the onset of extinction, since the 
heat transfer away from the flame is countered 
by the heat generated inside. This new obser- 
vation, therefore, suggests that in the setting of 
a “turbulent” flame in which the transitional 
effects are not dominant, and where the chem- 
ical reaction is describable by an Arrhenius 
model, the heat release results in a higher 
reactant conversion, even though the rate of 
mixing may be somewhat less. 

Based on our observation, we recommend 
that the effects of exothermicity be assessed by 
means of laboratory measurements. These 
measurements must involve a reacting system 
whereby the rate of reaction conversion is tem- 
perature dependent (unlike that employed in 
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Fig. 8. Scatter plot of product mass fraction versus the mixture fractions for Flame E and its comparison with the 
laminar diffusion flamelet model, x = 10 D. 


Ref. 39, and in which the large-scale mixing 
intensity is not significantly affected by the 
heat release (see Refs. 11, 12, 17-19, 25 for 
discussions on the effects of heat release in 
low speed combustion, and Refs. 14 and 20 for 
discussions of such effects in high-speed com- 
bustion). It is important to note that the issue 
would not be resolved by a mere comparison of 
the mixing characteristics between a reacting 
and a nonreacting layer. An appropriate analy- 
sis requires the consideration of two flows with 
the same hydrodynamic characteristics but with 
different magnitudes of the enthalpy of com- 
bustion. The extent of validity of our conclu- 
sion under more complex chemistry models 
can be determined by these experiments. 

CONCLUDING REMARKS 

A two-four compact finite difference algo- 
rithm is employed for direct numerical simula- 


tions of a spatially developing planar jet flame 
under nonequilibrium chemical conditions. 
With the high numerical accuracy provided by 
this algorithm, it is possible to perform direct 
simulations under physically realistic condi- 
tions of variable density and exothermic reac- 
tion. Here, the kinetic mechanism is assumed 
to obey the idealized model of the type F + 
O -» Products. 

The data produced by DNS are statistically 
sampled to assess the compositional structure 
of the flame. The results of this analysis are 
shown to be consistent with those of laboratory 
measurements, even with the assumption of an 
idealized chemistry model. It is shown that as 
the intensity of mixing increases, the magni- 
tudes of the ensemble mean and variance of 
the product mass fraction decrease, whereas 
those of the reactants increase. This conclu- 
sion is valid if the physical mechanisms by 
which the growth of instability waves is sup- 
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Fig. 9. Scatter plot of product mass fraction versus the reduced Damkohler number and its comparison with the 
laminar diffusion flamelet model. Mixture fraction in the range 0.45 < £ < 0.55, and * = 10£>. 


pressed are masked. This is implemented here 
by means of imposing a low amplitude random 
forcing at the entrance of the jet. The marginal 
and the joint pdfs of the reactants mass frac- 
tions, extracted from DNS data, show features 
in accord with laboratory data. That is, as the 
relative intensity of mixing is increased and the 
value of the Damkohler number is decreased, 
the pdfs show a lesser probability of product 
formation. The scatter plots of the instanta- 


neous product formation indicate that the lam- 
inar diffusion flamelet model is capable of 
predicting the compositional structure of the 
flame when the magnitude of the Damkohler 
number is sufficiently large. Under this condi- 
tion, the nonequilibrium effects are parameter- 
ized reasonably well by the instantaneous scalar 
dissipation rate in a simple laminar configura- 
tion with the same chemical characteristics. 
However, the flamelet library constructed in 
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Fig. 10. Streamwise variation of the total product mass fraction. 


this way does not yield satisfactory results at 
low Damkohler numbers. The results further 
indicate that in the setting of a “turbulent” 
flame, exothermicity results in an enhanced 
reactant conversion and an increased rate of 
product formation. This trend is not in accor- 
dance with previous DNS and laboratory 
results. The factors resulting in this discrep- 
ancy are (1) imposition of external forcing, and 
(2) the implementation of an Arrhenius kinetic 
model. Both of these factors are very impor- 
tant in diffusion flames in which the flow is 
fully turbulent and the chemical reaction rate 
can be assumed to obey the Arrhenius law. 
Therefore, we propose that in typical hydrocar- 
bon turbulent diffusion flames, the effect of 
heat release is to increase the rate of product 
formation even though the rate of mixing may 
be somewhat reduced. 
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Abstract— A family of Probability Density Functions (PDFs) generated by Johnson-Edgeworth Translation 
(JET) is used for statistical modeling of the mixing of an initially binary scalar in isotropic turbulence. 
The frequencies obtained by this translation are shown to satisfy some of the characteristics of the PDFs 
generated by the Amplitude Mapping Closure (AMC) (Kraichnan, 1989; Chen et aL, 1989). In fact, the 
solution obtained by one of the members of this family is shown to be identical to that developed by the AMC 
(Pope, 1991). Due to this similarity and due to the demonstrated capabilities of the AMC, a justification is 
provided for the use of other members of JET frequencies for the modeling of the binary mixing problem. 
This similarity also furnishes the reasoning for the applicability of the Pearson Family (PF) of frequencies 
for modeling of the same phenomena. The mathematical requirements associated with the applications of 
JET in the modeling of the binary mixing problem are provided, and all the results are compared with 
data generated by Direct Numerical Simulations (DNS). These comparisons indicate that the Logit-Normal 
frequency portrays some subtle features of the mixing problem better than the other closures. However, none 
of the models considered (JET, AMC, and PF) are capable of predicting the evolution of the conditional 
expected dissipation and/or the conditional expected diffusion of the scalar field in accordance with DNS. It 
is demonstrated that this is due to the incapability of the models to account for the variations of the scalar 
bounds as the mixing proceeds. A remedy is suggested for overcoming this problem which can be useful in 
probability modeling of turbulent mixing, especially when accompanied by chemical reactions. While in the 
context of a single-point description the evolution of the scalar bounds cannot be predicted, the qualitative 
analytical-computational results portray a physically plausible behavior. 


1 INTRODUCTION 

The problem of binary mixing in turbulent flows has been the subject of widespread 
investigations over the past two decades (Dopazo, 1973; Pope, 1979; Pope, 1985; Pope, 
1990; Givi, 1989; Kollmann, 1990). This problem has been particularly useful in assessing 
the extent of validity of the closures developed within this period for modeling of turbulent 
mixing by Probability Density Function (PDF) methods (Dopazo, 1973; Pope, 1976; Pope, 
1979; Janicka et aL, 1979; Pope, 1982; Kosaly and Givi, 1987; Norris and Pope, 1991). 
Usually the problem is considered in the setting of a spatially homogeneous turbulent flow 
in which the temporal evolution of the PDF is considered. In this setting, development 
of a closure which can accurately predict the evolution of the PDF has been the main 
objective of these investigations (for recent reviews see Pope (1990); Kollmann (1990); 
Givi (1989)). 

Computational experiments based on Direct Numerical Simulations (DNS) have proven 
very useful in evaluating the performance of new closures (Givi, 1989; Pope, 1990). The 
binary mixing problem is well-suited for DNS investigation, and current computational 
capabilities allow consideration of flows at sufficiently large Reynolds numbers in which 
the behavior of the models can be assessed (Eswaran and Pope, 1988; Givi and McMurtry, 
1988; McMurtry and Givi, 1989; Madnia and Givi, 1992). The results of all the previous 
work on DNS of the binary mixing problem portray a clear picture of the PDF evolution, 
at least at the single-point level. A successful closure is one which can predict all the 
stages of mixing, as depicted by DNS, from an initially binary state (total segregated) to 
a final mixed condition. 
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, A wT,f the mode Is developed in the literature, the recent Amplitude Mapping Closure 
(AMC) (Kra.chnan, 1989; Chen etal., 1989; Pope, 1991) has proven effective in producing 
a physically correct PDF evolution. In the application of this model to the problem of 
bmary mixing it has been demonstrated that the closure is capable of approximating a 
reasonably correct evolution at all stages of mixing (Pope, 1991). Namely, the evolution 

fr °. m , a " T t,al d ° ub e de ta PDF to an as y m P { °t'c Gaussian distribution. This is a trend 
which has been observed in DNS (Eswaran and Pope, 1988; Givi and McMurtry, 1988; 
McMurtry and Givi, 1989; Madnia and Givi, 1992) and also corroborated by experimental 
investigations (Miyawaki el al. , 1974). However, it is shown by Gao (1991)- O’Brien and 
Jiang (1991) that the PDF adopts an asymptotic Gaussian-Iike distribution ’only near the 
mean scalar value and the conditional expected dissipation does not correspond to that 
ot a Gaussian field everywhere within the composition domain. 

Our first objective in this work is to present another means by which the AMC can 
be viewed. It is demonstrated that in the binary mixing problem, this closure can be 
considered as a member of the family of frequencies generated by the method of Johnson- 
Edgeworth Translation (JET). In fact, it is shown that the result produced by the AMC is 
identical to that generated by one of the members of this translation. With this observation 
a justification for employing other simpler “assumed” frequencies is provided. Our second 
objective is to make a detailed examination of the conditional expected dissipation and 
me conditional expected diffusion of the scalar variable as predicted by the closures. 
This examination provides an effective means of demonstrating the deficiencies of these 
models in reproducing the correct physical behavior as depicted by DNS results. With 
the development of analytical relations for some of these closures, a remedy is suggested 
for overcoming the model deficiencies. 66 

/. 1 Outline 

In the next section, the problem of binary mixing and its solution via the AMC is briefly 
reviewed. In Section 3, the Johnson-Edgeworth Translation is introduced with a highlight 
on the mathematical constraints associated with its application for the modeling of the 
mixing problem. Due to the previously established similarity of the JET frequencies 
with those based on the Pearson Family (PF), the Beta density of the first kind is 
also Presented in this section. The PDFs generated by these three models (AMC, 
JET and PF) are compared against each other and also with data generated by Direct 
Numerical Simulations (DNS) in Section 4. The results for the conditional expected 
scalar dissipation, and the conditional expected scalar diffusion for all the closures are 
discussed, respectively, in Section 5 and in Section 6. In Section 7, some theoretical 
remarks pertaining to the evolution of the scalar in an isotropic field are presented With 
this presentation, the problems associated with all three closures become more clear, 
in Sections 2-7, the discussions are limited to those associated with the transport of a 
passive scalar from an initially symmetric binary state in isotopic turbulence. Therefore 
in Section 8 some discussions are presented of the applications of the models for treating 
more general problems. This paper is drawn to a conclusion in Section 9. 

2 BINARY MIXING PROBLEM 

We consider the mixing of a scalar variable cf> = <p(x,t) (x is the position vector, and t 
denotes time) from an initially symmetric binary state within the bounds <h < 6 < 6. In 
this section, we assume that the lower and the upper bounds of the scalar range"remain 

fixed (i.e. <t> u ,<f> ( are constant). Within this domain, the single-point PDF of the variable 
<fi at initial time is given by 

t = 0) = - 4> ( ) + 6(d> - 0„)J, 


(1) 
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15.0 



which obviously infers the following relations for the mean and the variance 


<<£>= (<f) u + 4>t)/ 2, < 4>' 2 >= o 1 = (2) 

Here, <> indicates the probability mean, a 2 denotes the variance, and the prime 
represents the instantaneous deviation from the mean. In isotropic incompressible 
turbulence, the evolution of the PDF is governed by the transport equation 


dPy d 2 {eP x ) 
dt d4 > 2 


4>t < 4> < 4> u 


(3) 


where e represents the expected value of the scalar dissipation with diffusion coefficient 
r, £(= r V<p ■ v<t>), conditioned on the value of the scalar 0(x, t), 

e = =< £\ (4) 

Equation (3) can alternatively be expressed by 


dP\ d(DPj) 
dt d(f> 


(f>f < <f> < (f> u . 


where D denotes the conditional expected value of the scalar diffusion 


(5) 


D - D(<p,t) =< rv 2 ^(x,r)>. 


( 6 ) 
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FIGURE lb The comparison of the PDFs as predicted by the models with DNS data, (b) a 2 = 0.079. 

The entire problem in determining the PDF, P u is associated with the unknown 
conditional expected dissipation, e, and/or the conditional expected diffusion D These 
two are related through Eqs. (l)-(6), ’ ' 56 


D{4>,t) = 


1 djePj) 

P\{<t>,t) d<j> - 


( 7 ) 


'T' "!. i,her ,he condi,ional mean dissipation nor the conditional mean 
diffusion are known (neither are their unconditional mean values). Their specifications 
require external information. ^ ,ons 

aPpli f ati °!l ° f ^ C ’ this external information is obtained in an implicit 
e ? plamed ' n det£ul b y P °P e ( 1991 X the AMC involves a mapping of the 
random field of interest 4> to a stationary Gaussian reference field <b, via a transformation 

■ \*a 1 ’ ? nCC r thlS relatlon is established, the PDF of the random variable 6, P x (6) 

is related to that of a Gaussian distribution. In a domain with fixed upper and lower 

Pnn!» noon the : corresponding form of the mapping function is obtained by 

pe (1991). The solution for a symmetric field with zero mean, < d> >= 0 d> = —6> 
is represented in terms of an unknown time t > w* 


x{4>0, T) - (f) u erf ^ ) ’ G( t ) = V cx p(2r) — 1. 

With this transformation, the PDF is determined from the physical requirement 

P\ (x(<k, r ), r)d x = P G (<t>o)d<f > 0 , 


( 8 ) 


( 9 ) 
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FIGURE lc The comparison of the PDF’s as predicted by the models with DNS data, (c) cr 2 = 0.00247. 


where Pc denotes the PDF of a standardized Gaussian distribution, i.e. Pc(<fa) = 
exp(— <Pq/2). A combination of Eq. (8) and Eq. (9) yields 


t) 



( 10 ) 


In these equations, the relation between r and the physical time t is unknown in 
the context of a single-point description. This relation can be obtained only through 
knowledge of the higher order statistical properties of the scalar field. For example, it 
is shown by Madnia et aL (1992); Frankel et ai (1992a) that the mapping closure yields 
the algebraic relation for the normalized variance, 


< <7 2 > (r) 
<< 7 2 > ( 0 ) 


2 

= — arctan 

7T 


( I 

\GVG 2 + 2 



( 11 ) 


(3), 

da 

where. 

a d 7 


(12) 


t>4>U m^U 

(0= / PM,t)e{<t>,t)d<t> =- (f>Pi (<t>,t)D(4>,t)d <f>. 

* <j>t J 4>t 


( 13 ) 



26 


R S MIIl.FR, S. II FRANK! 1., C. K MADNIA AND P GIVI 


3 JOHNSON-EDGEWORTH TRANSLATION 

h h p AM noo^rS S °T ° f ‘ hC b3SiC featUreS ° f the binar y mixin 8 P rob »™ as described 
by Pope (1991). Namely the inverse diffusion of the PDF in the composition domain 

from a double delta distribution to an asymptotic Gaussian distribution centered around 
< (J >, as a - 0 (or G -» oc). This asymptotic Gaussian distribution near the mean 
scalar value cannot be realized in any of the previous mixing models based on the so 

v eSCe a C ^ DiSf ; e nc^ n 1 C/D) Cl ° SUreS (Curl ’ I963; Janicka * < 1979; Pope, 
1982, Kosaly and Givi, 1987). And those modified C/D models which do yield such an 

^o^rigR^f Th'^'h P fi° Pe (198 ?l d °r> f S t PrediC * thC initia ‘ Stages 0f mbcin g correctl y 

(Kosaly, 1986). This deficiency of the C/D models in yielding asymptotic Gaussianity has 

AMC a (Pope V 199f) faCt ° r ^ reCCnt investigations resu,tin 8 in the development of the 

In the spirit of “mapping” to a specified reference field, it is speculated that there are 
perhaps other means of “driving” the PDF toward Gaussianity in a physically acceptable 
manner. In fact, this subject has been of major interest to statisticians and biometricians 
within the last century since the early work of Edgeworth (1907). The scheme was 
referred to by Edgeworth as the Method of Translation, and was later detailed bv 
Johnson (1949a). In today’s literature of statistics and biometrics, the scheme is known 

as Johnson-Edgeworth frequency generation, and has many applications in statistical 
analysis. 

t /he essent^ element of Johnson-Edgeworth Translation (JET) is similar to that of 
the AMC. Namely, it involves the transformation of the random physical field, here 6 
the^orm^ St3ndard Gaussian refer ence field by means of a translation (or mapping) of 

0o 


0 = 2 


1 . 7(0 J 


y(r _ 0) — 0 < 7 (0 < 7 (t —* oc) 


oo. 


(14) 


In this equation, the function 7 (t ) plays a role similar to that of G in the AMC. With 
an appropriate form for the function Z, the scalar PDF is determined from Eq. (9). 

or application in the problem of mixing from an initially symmetric binary state of zero 
mean within a fixed domain 4>t = -<j> u < 4> < 4> u , the appropriate JET must satisfy the 
following physical constraints: 3 


(i)Liin M Z(^)aff(^) 

7 


(iT‘)Lim (7 _ 00) Z (— ) « C0 O + 0(4) + ... 

(*>■) 2(— ) is an odd function with respect to the scalar mean for any value of a~. 

(iv) Z(^) is bounded and is a non-decreasing function of fo, and — d>„ < Z < 6 at 
all times. ~ - 

( 15 ) 

In these relations, H denotes the Heaviside function, and C is constant. Constraint (i) 
implies an initially symmetric and segregated binary state. The second constraint ensures 
an asymptotic Gaussian distribution for P\(4>) near the mean scalar value. Condition (iii) 
preserves the symmetry of the PDF around the mean value at all times, and constraint 
(iv) implies the boundedness of the scalar field, i.e. -0„ < 0 < 0 M . A function Z 



PROBABILITY MODELING IN TURBULENT FLOWS 


27 


which satisfies all the above constraints, is therefore expected to provide an acceptable 
means of approximating the PDF. An example is the Logit-Normal (or tanh -1 -Normal) 
distribution, as originally proposed by Johnson (1949a). For the symmetric problem 
within a fixed domain, this distribution is produced by a mapping of the form 1 

Z = <p u tanh . ( 16 ) 

With this mapping, together with Eq. (9), the PDF of the scalar adopts the form, 


P\(<t>J) 






tanh '(j - ) 

<Pu 


(17) 


It is easily verified that this frequency satisfies the physical constraints of the symmetric 
binary mixing problem (Eq. (15)). At t = 0, the PDF is approximately composed of two 
delta functions at 0 = ±4> u , and as t — ► oo the PDF adopts an approximate Gaussian- 
like distribution centered around the zero mean. These features are similar to those 
portrayed by the PDF generated by the AMC (Eq. (10)). 

This example demonstrates that with the satisfaction of the above indicated constraints, 
several other frequencies can be generated for effective modeling of the binary mixing 
problem. In fact, it is easy to show that the solution generated by the AMC can also be 

viewed as a member of the JET family. This is demonstrated by considering a translation 
of the form 

Z = *- erf (7) • (18) 

From Eq. (9), this translation yields the PDF 




57k exp r ( ^ _1) 


erf 




(19) 


This frequency can be termed the erf" '-Normal distribution and is identical to the form 
presented by Eq. (10). The difference is due to the terms containing G and 7. But this is 
unimportant since in the context of single-point statistics neither of the two parameters 
can be determined by the PDF. Therefore, with G = both expressions are equivalent. 

With this equivalence, the closed form relation for the variance of the erf" '-Normal 
distribution has the same algebraic form given by Eq. (11). It is easy to show that many 
other distributions can be generated to display similar characteristics. In the discussions 
to follow, we only consider the Logit-Normal and the erf" '-Normal distributions, the 
latter being identical to the distribution generated by AMC. 


Pearson Family, 

The similarity of the AMC and JET in generating equivalent PDF’s is also useful in 
explaining the applicability of the frequencies generated by the Pearson family (Pearson, 
895). For a bimodal distribution, a physically acceptable frequency is the Pearson 


I 


In recent literature, the Logit-Normal is usually expressed by the mapping Z = <j> u |2[1 + exp(0 o /7 )] -1 - l 
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FIGURE 2 Temporal evolution of the Logit-Normal PDE 

Type I, known as the general form of the “Beta density of the first kind”. This density 
is typically expressed in a fixed domain within the range 0 < <j> < 1, 

(20) 

Here B denotes the Beta function, and the parameters fi\ and are determined from 
the knowledge of the mean and the variance of the random variable. In a symmetric 

field within [0,1], < 4> >= 0\ — fh = 0, and thus the PDF is characterized by the 

variance alone. 

The similarity of the Pearson distributions and the JET frequencies is well recognized 
in the statistics literature (see Johnson (1949a)). Therefore, with the equivalence of the 
AMC and the JET as demonstrated above, it is not surprising that the Beta density 
and the AMC are also similar. This similarity, without a mathematical proof, has been 
recognized in previous works (Madnia et ai, 1991; Madnia and Givi, 1992). 

4 COMPARATIVE ASSESSMENTS 

The probability distributions obtained from the three frequency generation methods 
described above are all capable of providing a reasonable stochastic approximation of 
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FIGURE 3 a Temporal evolution of the centralized moments of the scalar variable as predicted by the 
models and the comparison with DNS data, (a) /14 vs. t*. 


the mixing problem from an initially binary state. Namely, an approximate double delta 
distribution at / =0, and an approximate Gaussian-like distribution as / — ► oo. The 
former can be realized in the limit of unity normalized variance <r 2 (/)/er 2 (0) = 1 . In a 
fixed composition domain, this corresponds to G = 0 for AMC (7 = 0 in JET), and 
to P = 0 for the Beta density. The latter is realized in the limit G,j,(3 —* 00 . The 
limiting Gaussian distribution for the AMC has been asserted t>y Pope (1991). For the 
JET, the criterion (ii) in Eq. (IS) guarantees this condition. For the Beta density, the 
assertion of an asymptotic Gaussian distribution in the limit of zero variance is established 
in elementary texts on statistics (e.g. Casella and Berger (1990)). At the intermediate 
stages, however, the PDF’s are not identical. It is easily verified by Eqs. (10), (20) that 
the AMC and the Beta distributions become constant (P|(<£) = constant = \4> u ) for 
G = (3 = 1. However, the Logit-Normal PDF does not yield a uniform distribution at 
any stage of its evolution. Also, as indicated by Johnson (1949a) it is not possible to 
provide a closed form algebraic expression similar to Eq. (11) for the variance of the 
Logit-Normal distribution. 

In order to make comparative assessments of the models, the frequencies generated 
by the three methods (AMC, JET, and PF) are compared with each other, and also with 
PDF’s generated by Direct Numerical Simulations (DNS). The DNS procedure is similar 
to that of previous simulations of this type. Since these simulations are not the major 
focus of this paper, only a brief outline of the procedure is described; for a detailed 
discussion we refer the reader to Madnia and Givi (1992). The subject of the DNS is 
a three-dimensional periodic homogeneous box flow carrying a passive scalar variable. 
The initial scalar field is composed of square waves with maximum and minimum values 
of 1 and 0, respectively. These limiting values are arbitrary, and can be translated to 
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FIGURE 3b Temporal evolution of the centralized moments of the scalar variable as predicted by the 
models and the comparison with DNS data, (b) /r 6 vs. t*. 


a PPropri a te 4 > u , for comparison to each of the models. At time zero, in most regions 
of the flow, the scalar adopts these limiting values (with equal probability), and the 
spatial regions between the initial maxima and minima are smoothed by a polynomial fit. 
The degree of this smoothing is minimized, but limited by the computational resolution, 
to ensure an approximate initial double delta distribution. The computational scheme 
is based on a spectral-collocation procedure using Fourier basis functions developed 
by Erlebacher et al. (1990a); Erlebacher et al. (1987); Erlebacher et al. (1990b). The 
hydrodynamic field is assumed isotropic and is initialized in a similar manner to that by 
Erlebacher et al. (1990a); Passot and Pouquet (1987). The code is capable of simulating 
flows with different levels of compressibility (Hussaini et al., 1990). Here, only the 
results obtained for a low compressible case are discussed. The resolution consists of 
96 collocation points in each direction. Therefore, at each time step 96^ is the sample 
size for statistical analysis. With this resolution, simulations with a Reynolds number 
(based on the Taylor microscale) of Re\ as 41 are attainable. The value of the molecular 
Schmidt number is set equal to unity. 

As indicated in Section 1, in order to compare the model predictions with DNS results 
a matching is required of the higher order statistics of the field as generated by each 
method. Here, this matching is done through the variance of the conserved scalar. These 

results are presented in Fig. 1. This figure indicates that at initial times, as 1, all 

the PDF’s are approximately composed of two delta functions at <p = 0,1 indicating the 
initial binary state. At longer times, the PDF’s evolve through an inverse-like diffusion 
in the composition space. The heights of the delta functions decrease and the PDF’s 
are redistributed at other <f> values within the range [0,1]. At very long times, the PDF’s 
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FIGURE 4a Comparisons of the normalized conditional dissipation as predicted by the three models with 
the DNS data, (a) <r 2 = 0.079. 

become asymptotically concentrated around the mean value in a manner that can be 
approximated by a Gaussian distribution. 

An interesting feature captured in Fig. 1(b), is the capability of the Logit-Normal 
distribution in depicting a subtle behavior in the frequency distribution. This feature is 
the double “hump” characteristic of the DNS data at intermediate times and cannot be 
realized by the AMC or the PF generated frequencies. All the previous DNS results 
including those of Eswaran and Pope (1988); Givi and McMurtry (1988); Pope (1991) 
portray this feature. The PDFs generated by the AMC, and the Beta distributions 
adopt a constant value (of 1/2) when a 2 = ^ (for 0 < (f> < 1). This corresponds to 

G = 1,7 = a/ 2, = 1. This uniform distribution is not exactly realized in any previous 

or present DNS results. Therefore, it can be speculated that in the absence of a better 
alternative, the Logit-Normal distribution may provide the simplest means of providing an 
assumed distribution for the statistical modeling of the symmetric binary mixing problem. 
The complete evolution of the Logit-Normal PDF is shown in Fig. 2. 

Further quantification of the agreements noted above are made by comparing the higher 
moments of the scalar field. This comparison is made in Fig. 3. In this figure, results are 
presented for the temporal variations of the kurtosis (m) and the superskewness (/i 6 ) 
of the scalar variable 4>. For the Beta density, the higher order moments are obtained 
analytically based on the knowledge of the variance. For the AMC, the analytical- 
numerical results by Jiang et al. (1992) are used, while for the Logit-Normal PDF the 
moments are calculated strictly by numerical means. This figure shows that initially, all 
these moments are close to unity, and monotonically increase as mixing proceeds. For all 
the models, the magnitude of the moments asymptotically approach the limiting values 
of 3 and 15, respectively, corresponding to those of a Gaussian distribution. The DNS 
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results are in good agreement with the model predictions at all times. However, due to 
obvious numerical difficulties, the simulations could not be continued until the variance 
approaches zero identically. 


5 SCALAR DISSIPATION 

The results presented above indicate a good agreement between the model predicted 
single-point statistics (PDFs and high-order moments) and the DNS data at all the 
stages of mixing. These results also suggest an approximate asymptotic Gaussian state 
for all the closure PDFs and those of the DNS. Here, it will be demonstrated that the 
agreement between the DNS and the model predictions is very good at the initial and 
the intermediate stages of mixing. However, the agreement worsens at the final stages 
Also it will be shown that none of the closures yield “exact” Gaussian distributions at the 
final stages of mixing. In doing so, it is useful to note that a Gaussian PDF is defined, 
and is only valid, for an unbounded domain. The frequencies generated here, are all 
defined within & fixed and finite domain. For AMC, it has been established (Gao, 1991- 
O’Bnen and Jiang, 1991) that the finite boundary size at the initial time “maintains” 
its influence at all the subsequent stages of mixing. In other words, the PDF adopts a 
Gaussian distribution in the limit of zero variance only near the mean value of the scalar. 
In order to show the departure from Gaussianity at scalar values away from the mean 
the conditional expected dissipation of the scalar field is considered. 

Given the PDF, as is the case here, Eq. (3) can be used to determine the expected 
conditional dissipation. It has been shown by Girimaji (1992) (and will be discussed in 
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FIGURE 5a Temporal evolution of the conditional expected diffusion normalized with the total dissipation, 
(a) erf - 1 -Normal. 


detail in Section 7) that for a valid PDF within a defined range, -<f> u < <f> < (j> u , the 
expected conditional scalar dissipation is given by 

(21) 

where F denotes the cumulative distribution function (CDF) 

0= / (22) 

J-+U 

With Eq. (21), the expected conditional dissipation can be evaluated for a given PDF. For 
example, for a Gaussian distribution of zero mean, P g (4>,<t 2 ) = exp(— ^), — oo = 

-<f> u < <f> < <t> u = oo, with a non-stationary variance, a 2 = <r 2 (t), it is easily shown that, 

, — oo < <f> < oo 

(23) 

Noting that ^ is an independent variable (of /), and evaluating the derivatives on the 
RHS of Eq. (23) yields, after some simple manipulations, 

/ , \ da 

e(<M) = constant = — < 7 — , 


e(<M) = 
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L Pg(4>, o -2 ) 


1 ft 

dt \2 
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)+ 


\/2i r 
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Temporal evolution of the conditional expected diffusion normalized with the total dissipation 
(b) LMSE. r 


with the implications (derived from Eqs. (12)-(13)), 


£(<M) , . 

— — - = constant = 

«(0 


1, 


(25) 


at all times for all <f> values in the range -oc to +oo. Equations (24-25) indicate the 
independence of the conditional scalar dissipation and the composition domain for a 
Gaussian field. These results have also been obtained by Gao (1991); O’Brien and Jiang 
(1991) by following a different mathematical procedure. 

The conditional expected dissipation predicted by the models can be obtained by 
following a similar course. For the AMC and the PF distribution, the conditional 
dissipation fields have been obtained by Gao (1991); O’Brien and Jiang (1991) and by 
Ginmaji (1992), respectively. For the purpose of the discussions to follow, these results 
are presented here in a different form for all three closures. For the erf - '-Normal 
distribution, the instantaneous CDF is given by 




-K 


1 + erf 




) 


(26) 


Therefore, with Eq. (21), the conditional dissipation can be expressed in terms of the 
corresponding PDF, 


f(<M) = ~ 


1 d 
PM,t)dt 


tic 


erf 




d<fr + j (0 + ^ 




(27) 
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Again, with an independence of </> and t this equation reduces to. 


e(4>' 0 = 


1 OH 
P l (4>.t)~dT' 


(28) 


where, 

4> u r rf ~' { £ ) 7 , 

H=-= erf(-j=z)exp(-z 2 )dz. (29) 

V n J —<t> u v2 

For a PDF within a fixed domain, the integration procedure becomes simplified by 
evaluating the time derivatives inside the integral. In this way, the results can be expressed 
analytically. After some manipulations, 


dH 

~dt 


— y/24>u jj- 
?t(2 + 7 2 ) 


exp 


{ 


7 

-< , + 7 > 


erf ~'(T") 

Vu 


(30) 


and, therefore 




tt7(2 + 7 2 ) 


exp 



<Pu 



(31) 


From this equation, the total dissipation is obtained by direct integration of the conditional 
mean dissipation field. The results, after significant algebraic manipulations yield 


€(/) = 



7r(2 + 7 2 )\/4 + 7 2 


(32) 


g(0.O 

<0 


,N 


1 + sin 



erf-'(-f) 


(33) 


In the form presented above, Eqs. (31)-(33) portray several insightful features of the 
solution. First, Eq. (33) indicates that the conditional dissipation is always dependent 
on the magnitude of the scalar, and it maintains the same self-similar functional form of 


dependence exp 


{ 


-2 



This has been previously indicated by Gao (1991); 


O’Brien and Jiang (1991). Here, the amplitude e(<p = 0, / ) can be conveniently expressed 
in terms of the variance decay, which is very useful for ftirther manipulations. Second, 
it is interesting to note the similarity of Eqs. (31) and (33) with the results obtained for 
the instantaneous dissipation of Fickian mixing of a conserved scalar in laminar non- 
homogeneous flows (such as the typical shear flows (Spalding, 1961; Linan, 1974; Peters, 
1984)). This similarity further asserts the “permanent” influence of the boundaries 
since in non-homogeneous mixing, the scalar bounds are “fixed” due to the physical 
constraints. Finally, Eq. (32) suggests an infinitely large dissipation at time zero, i.e. 
when cr 2 (r)/cr 2 (0) == 1, and the asymptotic behavior 
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FIGURE 6a Comparisons of the normalized conditional diffusion as predicted by the three models and 
the LMSE closure, with the DNS data, (a) a 2 = 0.079. 


This limiting behavior near zero indicates the Gaussianity of the PDF only at the mean 
value of the scalar. 

Following a similar procedure, the conditional expected dissipation can be obtained 
for the other closures. For the Beta density in the range 0 < <f> < 1, the final results can 
be expressed as 

'**»**)• us) 

where, I denotes the Incomplete Beta Function (Abramowitz and Stegun, 1972). For 
the Logit-Normal distribution, the corresponding form is 


e(<M) = 


1 


d 


Pi(<j>,t)dt 


{1jC-[*‘G3)M 


(36) 


Neither of the equations (35-36) can be simplified further. Therefore, in order to evaluate 
the conditional expected dissipation (and the total dissipation), these equations must be 
evaluated numerically. 

In Fig. 4, the evolution of the conditional expected dissipation (normalized by the 
total dissipation) is presented for the models and the DNS data. This figure shows the 
similarity of the conditional expected dissipation for all of the models. The bell shape 
distribution is evident in all the figures with a maximum amplitude near the mean value. 
Also, as the variance decreases and the PDF becomes concentrated near the mean, the 
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FIGURE 6b Comparisons of the normalized conditional diffusion as predicted by the three models and 
the LMSE closure, with the DNS data, (b) a 1 = 0.013. 

amplitude tends to unity. This shape is typical of that observed in previous DNS results 
of Eswaran and Pope (1988); Nomura and Elgobashi (1992). 

The results in Fig. 4 show the <f> dependency of the results at all the stages of mixing. 
That is, the PDF asymptotically adopts an apparent Gaussian-like distribution only near 
the mean value of the scalar, and the conditional dissipation does not become independent 
of the scalar everywhere. For the AMC, this has been discussed by Gao (1991); O’Brien 
and Jiang (1991). Considering the similarity of the three models, it is therefore concluded 
that all three models yield the same characteristics. These results also suggest a poor 
agreement between the model predicted conditional expected dissipations and the DNS 
data. Note that at the initial stages of mixing, the predicted results compare well with 
DNS data. However, with mixing progression, at smaller variance values, the agreement 
is only good near the mean scalar value and worsens near the bounds of the composition 
domain. This, as described above, is due to the permanent influence of the scalar 
boundaries at all the stages of mixing. 

6 SCALAR DIFFUSION 

Albeit directly related to the conditional expected dissipation (Eq. 7), it is useful to 
examine the behavior of the conditional expected diffusion in light of the discussions 
above. Given the PDF, again within the fixed range — <f> u < <f> < (f> u , the conditional 
expected diffusion can be determined from 


1 dF 




) dt ‘ 


(37) 
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This equation is very useful in illustrating the properties of the PDF For example, for a 
Gaussian distribution within an infinite domain 


dF _ 4> da 2 ft 2 

~ ~72^~dT CXp{ ~2^ ) 


(38) 


and consequently 


£& 0 = ~4> 

e(0 a 2 {t)' 


(39) 


It is noted that Eq. (39) is in accord with the Linear Mean Square Estimation (LMSE) 
closure (O’Brien, 1980). v 

The mean conditional diffusion can be determined for the three models considered. 
For the erf '-Normal PDF with zero mean 


OF 

dt 


1 d-y 

W,7J exp 



<Pu 
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(40) 


Again with explicit equations for the total dissipation and the variance, it is possible to 
obtain an algebraic expression for the conditional expected diffusion. The results after 
substantial algebraic manipulations yield 


<t) 




1 + sin 
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Hto i \ 

[2^(0)J ) 


1 — sin 
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' l£lLLl 


exp< - 


{■ 


<Pu 


erf '(r - )- ( 41 ) 


In this form Eq. (41) is very pleasing since it does reveal the (t ,</>) separability, and thus 
the self-similarity, of the diffusion field. The terms inside the parenthesis on the RHS are 
time dependent, whereas the remaining terms depend explicitly on 0 only. As indicated 
by O’Brien and Jiang (1991), this separability cannot be easily deduced from Eq. (5), 
but is possible with the analytical procedure followed above. The temporal evolution of 
the conditional expected diffusion for the erf-' -Normal distribution, and its comparison 
with that of the LMSE closure is presented in Fig. 5. 

By following the procedure above, analogous expressions are obtained for the other 
two closures. Namely, 


t ) _ 2 < fu dj_ 

c(0 7 da 2 


tanh-'(-^), 


(42) 


for the Logit-Normal distribution of zero mean, and 

WO = J_W) 

e(0 Pi (<M) da 2 (43) 

for the symmetric Beta density within [0,1]. Equations (42)-(43) cannot be simplified 
further due to the lack of an explicit analytical relation for the variance of the Logit- 
Normal distribution (Johnson, 1949a), and the unknown analytical form of the derivatives 
of the Incomplete Beta Function. 
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FIGURE 7 The temporal variations of <t> max and 4> mjn generated by DNS. 

In Fig. 6, results are presented of the conditional expected diffusion as predicted by the 
three models, and also that of the LMSE closure. In these figures, the DNS data are also 
provided at several variance values. The similarity of the modelled results are once again 
revealed in these figures, which is expected in view of the PDF similarities. At all times, 
the conditional diffusion field has an odd distribution near the mean scalar value. On 
the right half of the composition domain, all three closures yield a monotonic decrease 
of D to an instantaneous minimum, and then a monotonic increase to zero at the upper 
bound of the scalar. The location and the magnitude of the instantaneous maxima and 
minima is not the same for the three closures. Also, as Eqs. (41)-(43) indicate, the 
zeroes of D can only be realized at (f> = 0, ±<£ u . At the initial times, i.e. large variances, 
all three closures agree reasonably well with the DNS data. This agreement is better 
for the three models than for the LMSE closure. However, as the variance becomes 
smaller, the agreement between the model predictions and the DNS data worsens. It 
is noted that as the variance becomes small, all the closures yield a Gaussian-like PDF 
near the mean value of the scalar. This is shown in the figures near 0=<<^>(=i 
for DNS), where the predicted results are in accord with the LMSE closure, i.e. linear 
profiles of similar slopes. In this region, the results are also in accord with DNS data for 
all the closures including the LMSE. However, again, the predicted results deviate from 
the DNS data away from the mean value. It is clearly noted that the DNS generated D 
values do not go to zero at the scalar bounds. 

7 EVOLUTION OF THE SCALAR FIELD 

The problems described at the conclusion of the previous two sections stem from a 
lack of capability of all of the models in accounting for the variations of the scalar 
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FIGURE 8 The Beta density (Pearson Type II) for a domain with moving boundaries, and 0 = 0.1. 

bounds as the mixing proceeds. For all three models, the PDF is always defined within a 
fixed range through its course of evolution. It is easy to show that both the conditional 
expected dissipation and the conditional expected diffusion are correctly predicted by all 
the models near the mean scalar value. For the erf^Normal distribution, this is evident 
from Eqs. (33) and (41) and can be also shown by analyzing the behavior of Eq. (31) 
near the region <j> % 0, as the variance becomes small. Noting that 
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0 yfH 4> 
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<pu 2 <p u 

(44) 

and, from 

Eq. (11) 
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Following the same procedure, it is derived 

Fim ( , ->oo) = (47) 

Due to the similarity of the three closures, it is reasonable to expect similar behaviors 
for the other models as well. Equations (46)-(47) indicate a Gaussian-like distribution 
near the mean <f> v 0 (Eqs. (24,39)). This is in accord with the DNS data. However, at 
distances away from the mean value the predicted results do not correspond to that of 
a Gaussian field. Neither do these results agree well with DNS data. The deficiency of 
the models in predicting the DNS results is made clear by considering the bounds of the 
scalar field as the mixing proceeds. This is demonstrated in Fig. 7, showing the temporal 
decay/growth of the scalar maxima/minima obtained by DNS. This trend is consistent 
with physical intuition, but is not incorporated into any of the three models. In the 
AMC and the JET generated frequencies, due to the nature of the translation Z(<£ 0 , 0 
and the constraints imposed by Eq. (15), the scalar is always bounded within the same 
range. This problem is also encountered in the PF, in that Type I and II distribution 
families are always defined within the same fixed domain regardless of the magnitude of 
the variance. 

With the examination of the PDF transport equation, it is shown that the physics 
of the problem requires the migration of the scalar bounds toward the mean value 
as the mixing proceeds. That is, the instantaneous values of the scalar minima and 
maxima change with mixing progression. To demonstrate this, again consider a symmetric 
field with a PDF, P\(<j>,t), defined within the time-dependent domain of zero mean, 

< ^ f [ < / , min(0 = ~4>max(t ), 4>max(t )]• At all the stages of the evolution, the PDF must satisfy 
the physical requirements 

/*$max(* ) 
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The first of Eq. (48) requires 
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Evaluating this integral via Leibnitz’s rule, and making use of Eqs. (3)-(7), it is shown 
that 
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0max(O) — (f) (t , c/> mm (0) — (pf = — 0 W . (50) 

Following the same procedure for the second of Eq. (48), yields the obvious requirement 

f*0max (0 

D(<f>,t)P\(4>,t)d(t> = 0. (51) 


/ 

J<h, 




The third part of Eq. (48) yields Eq. (12), and 


0 £ (<Amax(0t 0 = ?\ (0min(O> O e (^min(f )* 0 = 0, 


P \ )■ 0 


‘/’min 

Jr 


^Wmi„(r),r) 


= 0. 


(52) 


The remaining parts of Eq. (48) yield higher order statistical information pertaining 
to the inner integrated evolution of the conditional expected dissipation and diffusion, 
and their relation with the higher central moments. With an additional assumption of a 
nonzero PDF within the region of its definition, that is by defining 0 max (/) and 
as the extreme locations with nonzero PDF, a combination of Eqs. (50) and (52) yields 


£ (^max(0; 0 — (O’O = 0, 


-Sf = ^Kax=D(0 max (r),r), 

^0min | , r> / , 

^ [ 0min = /), 

^max(O) = </»„, <Amin(0) = (53) 

Equation (53) indicates that with fixed boundaries, the conditional dissipation would 
adopt a zero slope at the boundaries and the conditional diffusion would also be zero 
there. Flowever, Fig. 7 indicates that in a physical situation the boundaries are not fixed 
and move inwards as the mixing proceeds. It is interesting to note that this problem 
is not observed in the numerical results obtained by the C/D type closures. That is, 
while the C/D closures are not capable of predicting the PDF evolution in accord with 
DNS data, they do have the mechanism for shrinking the bounds of the composition 
space. Obviously, in the context of single-point description without the knowledge of 
the dissipation field, it is not possible to determine a priori the temporal bounds of the 
scalar field. Therefore, the closures can be modified only by making further assumptions 
in describing this transport. For a general case, the JET frequencies can be generated 
by the original form proposed by Johnson (1949a) 


^ ) — A(r)z + p(Oi (54) 

where the additional parameters X(t ) and g(t ) provide the extra degrees of freedom in 
order to account for the variations of the instantaneous boundaries of the composition 
domain. For the PF, the problem can be overcome, for example, by considering a 
“four-parameter Beta distribution” 
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FIGURE 9a The comparison of the conditional expected scalar dissipation normalized with the total 
dissipation with DNS data as predicted by the AMC with the scalar bounds determined from the DNS results, 
(a) a 1 = 0.079. 
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with the extra two parameters being 0 m i n (O and </>max(0- F° r a symmetric PDF in the 
range [0,1] = [</> m i n (f = 0) = (j)e,<f> max (t = 0) = </>„]; therefore, the variance decay 
can be influenced by increasing /?, and/or by decreasing the scalar range A <j>(t) = 
</>max(0 - </>min(0' The former recovers the well-known two-parameter Beta distribution 
(Pearson Type II), while the latter is approximately equivalent to the LMSE closure 
(O’Brien, 1980). This latter case is presented in Fig. 8 showing a symmetric Beta density 
with P(t ) = fixed = 0.1. Note that as the mixing proceeds, the variance decays but the 
PDF preserves its initial approximate double delta shape. In a physical problem, the 
situation is somewhere between these two limiting cases. The exact situation depends 
on the characteristics of a particular mixing problem. 

The discussions above suggest that in order to predict the final stage of mixing correctly, 
the effects of mixing on the shrinkage of the domain must be taken into account. To 
demonstrate this point, the results shown in Fig. 7 can be incorporated into the mixing 
models to determine the evolution of the conditional expected dissipation and diffusion. 
This is done here only for the erf - '-Normal, and the results of the conditional expected 
dissipation are shown in Fig. 9. In the calculations resulting in this figure, analytical 
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FIGURE 9b The comparison of the conditional expected scalar dissipation normalized with the total 
dissipation with DNS data as predicted by the AMC with the scalar bounds determined from the DNS results, 
(b) <r 2 = 0.013. 

solutions are not possible for the moving boundary case. This is demonstrated by the 
equivalent form of Eq. (29) 


„ ( , ) = z¥£ b f^} dz + l 

K V2n Jo 1 + ^ 2 


( 4> + <t> max (0). (56) 


With this equation, therefore, the effect of the temporal variation of the PDF on the 
conditional dissipation is through the dH /dt term in Eq. (28). This term has the form 
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The first term on RHS of Eq. (57) is the same as that in Eq. (30), and the effects of moving 
boundaries manifest themselves through the second term. This term cannot be evaluated 
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analytically. However, Eqs. (56)-(57) show that due to the direct dependence on 
the conditional dissipation does not retain its original functional dependence, suggested 
by Eq. (33). Also, Fig. 9 shows that the effect of the moving boundaries is to force the 
conditional expected dissipation to zero at the current scalar maxima/minima. Therefore, 
the predicted results compare much better with DNS data than those presented in Fig. 
4. Due to the similarity of the PDF’s, it is expected that the other two closures would 
also behave in the similar fashion. 

The influence of boundary encroachment is also sensed in the conditional expected 
diffusion field. For the erf - -Normal scalar PDF, the equivalent form of Eq. (41) is 
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with an average dissipation of 
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Equations (58)-(59) show the influence of the boundary movement through the last term 
on the RHS of both these equations. With these additional terms, the normalized form 
similar to Eq. (41) is not very useful, and Eqs. (58)-(59) are evaluated numerically. 
The equivalent of Eq. (58) for the Logit-Normal and the Beta density are, respectively. 


D(<t>,t) 


— IT tanh ~‘ (— W) + 

^maxT dt \^m**(0/ ax(0/ dt 


(60) 


w °- 7 Sol('{.£Sfe} Wfl ) (6,) 

An interesting characteristic displayed by Eqs. (58) and (60) is the influence of the 
terms containing the temporal derivative of Note that at the boundaries, Le. 

<f> = 4>msu(t), the first term on the RHS of these equations vanishes, but the last term 
prohibits the conditional expected diffusion from going to zero. This is in accordance 
with the DNS data as shown in Fig. 6. In order to demonstrate this more clearly, results 
are presented in Fig. 10 of the conditional expected diffusion predicted by the erf -1 - 
Normal, with the input of the variance and the scalar bounds from DNS. A comparison 
between this figure and Fig. 6 show the influence of the boundary movement, and a 
better agreement between the model predictions and the DNS data. This agreement 
is more pronounced at scalar values away from the mean. Near the mean value, the 
influence of the boundary migration is slight, as also indicated by Eqs. (58) and (59). 


8 DISCUSSIONS AND APPLICATIONS 

In previous sections, a rather detailed discussion was presented of the problem of scalar 
mixing from an initially symmetric binary state. These discussions were primarily intended 
to provide a means of assessing the differences between the currently available tools in 
probability modeling of the scalar mixing problem. This problem is of significant interest, 
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FIGURE 10a The comparison of the conditional expected scalar diffusion normalized with the total 
dissipation with DNS data as predicted by the AMC and the LMSE closures with the scalar bounds 
determined from the DNS results, (a) a 2 = 0.079. 

considering the extent of previous works focused on its analysis (Pope, 1976; Pope, 1979; 
Pope, 1982; O’Brien, 1980; Dopazo, 1973; Kosaly and Givi, 1987; Pope, 1991; Gao, 
1991; O’Brien and Jiang, 1991; Nomura and Elgobashi, 1992). The results obtained here 
are particularly useful in highlighting some of the deficiencies of these closures, and in 
suggesting future research towards overcoming these drawbacks. There are, however, 
many other physical problems that are not subject to the restricting conditions imposed 
in these analyses. In this section, therefore, some discussions are presented as to the 
practical implications of these models, together with some speculations on their extensions 
for future applications. 

Perhaps one of the most important practical applications of the closures considered here 
is the treatment of reactive flow phenomena. In fact, the most important advantage of 
scalar PDF methods is due to their applicability in the modeling of turbulent combustion 
(Pope, 1979; Pope, 1985; Pope, 1990; Kollmann, 1990; O’Brien, 1980). The results 
generated here can be used directly in the modeling of mixing controlled homogeneous 
chemically reacting systems. Namely, in examining the compositional structure of a 
reacting system under chemical equilibrium, or in determining the limiting rate of 
reactant conversion in a simple chemistry of the prototype Fuel + Air — ♦ Products. 
The determination of this rate has been the subject of extensive investigations over 
the past forty years (see Hawthorne et al. (1949); Toor (1962); Williams (1985)). It is 
now well-established that in a mixing controlled binary irreversible reaction of this type, 
the statistics of the reacting fields can be related to those of an appropriately defined 
conserved scalar (such as <f>) (Bilger, 1980; Toor, 1975; Williams, 1985). Therefore, the 
frequencies generated herein can be utilized for estimating the statistics of the reacting 
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field with an infinitely fast chemistry model in a homogeneous flow with initially segregated 
reactants under stoichiometric conditions. Albeit very restricting, this problem is of great 
practical importance for modeling and design of batch mixers and plug flow reactors 
in which these conditions prevail (Toor, 1975; Brodkey, 1981; Dutta and Tarbell, 1989). 
Madnia et al (1991); Madnia et al. (1992) have shown that with the erf - '-Normal 
(AMC) and the Beta density models, this rate can be predicted by simple analytical 
means. For the Logit-Normal density, a complete analytical solution cannot be obtained 
and determination of the statistics requires numerical integration of the PDF. The results 
generated by these closures agree with DNS data better than those obtained by means 
of the C/D closures (McMurtry and Givi, 1989), or other models previously available 
in the chemical engineering literature (Dutta and Tarbell, 1989) (see Givi (1989) for a 
review). Also, the results provided by the AMC (Frankel et al ., 1992a) are shown to 
compare well with experimental data on plug flow reactors if the additional information 
pertaining to the evolution of the scalar length scale is accurately provided. 

The most obvious issues in regard to the applications of these models are associated 
with their extension for the treatment of (1) non-symmetric binary scalar mixing, (2) non- 
binary scalar mixing, (3) multiple scalar mixing, and (4) non-homogeneous mixing. The 
first problem constitutes a more general form of the binary mixing problem and is also 
important for the analysis of non-stoichiometric reacting systems. The second problem is 
appropriate for the analysis of other mixing systems in which the initial conditions are not 
of a two-feed configuration. The third problem is of interest in reacting systems in which 
the transport of a passive scalar (like <p ) is not sufficient for a complete analysis. For 
example, any reacting system under non-equilibrium conditions. Finally, the importance 
of the fourth problem is obvious in view of the fact that the flow within most practical 
mixing devices cannot be assumed homogeneous. 

In regard to the first issue, all of the three closures considered here can be used for the 
probability modeling of scalar mixing within a fixed scalar domain. The use of the AMC 
is straightforward, but the mathematical procedure is somewhat complex (Madnia et al., 
1992). The Pearson frequencies can be generated easily for non-symmetric problems. 
In this case, the Pearson Type I provides a reasonably accurate representation of the 
scalar field regardless of the degree of asymmetry of the PDF (Frankel et al, 1992b; 
Madnia et al, 1992). The use of the JET in this regard is most difficult, since closed 
form analytical expressions are not available for the variance of the scalar by which the 
PDF can be characterized (Johnson, 1949a). In treating these problems, therefore, the 
first two models can be more readily employed and subsequently used for the treatment 
of mixing controlled reacting systems under non-stoichiometric conditions. In fact, as 
demonstrated by Madnia et al. (1992) the solution of the non-symmetric form of the 
AMC and the Beta density provide a very good means of predicting the limiting rate 
of reactant conversion in homogeneous reacting flows. However, it should be indicated 
that with both models the problem associated with the scalar bounds still exists and must 
be dealt with as discussed in Section 7. 

In addressing the second issue, it is obvious that the AMC is more appropriate than 
the other closures for simulating the mixing problem from an initially “arbitrary” state. 
The extension of JET and PF for treating multi- (higher than bi-) modal distributions 
have been reported in statistics literature. However, as the degree of modality of the 
PDF increases the procedure becomes more complex and not suitable for practical 
applications. Fortunately, in most mixing problems in simple flows, i.e. homogeneous 
turbulence and turbulent shear flows, the PDF exhibits strong bimodal features (Madnia 
et al, 1992; Frankel et al, 1992b). In those cases, the use of the Beta density can be 
justified. In fact, in non-homogeneous flows it is easier to use this density, at least until 
further development of the AMC for practical applications (see Frankel et al. (1992b) - 
Gaffney et al. (1992)). 
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FIGURE 10b The comparison of the conditional expected scalar diffusion normalized with the total 
dissipation with DNS data as predicted by the AMC and the LMSE closures with the scalar bounds 
determined from the DNS results, (b) a 2 = 0.013. 


The extension of all of the three models in describing multi-scalar mixing is possible. 
The problem naturally falls within the realm of the multivariate statistical analyses. In these 
analyses, the implementation of the AMC is relatively straightforward since it provides 
a transport equation for the joint PDFs of the scalar variable in a multivariable sense 
(Pope, 1991). However, it is not presently clear how to devise an efficient computational 
procedure, typically based on Monte Carlo methods (Pope, 1981), for the numerical 
treatment of these equations. Some work in this regard is currently under way (Valino 
and Gao, 1991). The extension of assumed distributions based on the Beta density for 
treating multi-scalars is more straightforward but less trivial to justify. The most obvious 
means is to implement the multivariate form of the PF. The direct analog of the Beta 
density is the Dirichlet frequency (Johnson, 1987; Narumi, 1923; Johnson and Kotz, 
1972). The application of this density in modeling of multiple species reactions has been 
discussed by Girimaji (1991a); Girimaji (1991b); Gaffney et al. (1992). However, the 
use of the Dirichlet frequency cannot be justified for modeling of reacting flows in a 
general sense (Frankel, 1992). Finally the extension of the JET in generating multivariate 
frequencies has been reported in statistics literature since the subsequent work of Johnson 
(1949b). As one may suspect, the procedure is more complex, and the same reservations 
as those associated with the Dirichlet distributions apply. 

All of the models considered here can be extended for the analysis of non-homogeneous 
mixing (and reacting) systems. Obviously, in most cases, the problem requires numerical 
integration of the appropriate conservation equations. For instance, the AMC can be 
invoked in the mixing step of a fractional stepping procedure, similar to that of typical 
Monte Carlo procedures (Pope, 1981). The PF densities ( e.g . Beta or Dirichlet) and JET 
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generated frequencies require modelled transport equations for the first two moments 
and cross moments of the scalar field. These equations, “hopefully”, include the essential 
information pertaining to the spatial inhomogeneity of the flow. Naturally, the PDF is 
not generally symmetric, and must be determined from the knowledge of the parameters 
and the local ^ma*(0»^min(0 values. With this information, all the higher 
order statistics of the scalar field can be determined. In regard to this last issue, it must 
be indicated that the Beta density has been extensively used for the modeling of non- 
homogeneous reacting systems (e.g. Rhodes (1975); Jones and Priddin (1978); Lockwood 
and Moneib (1980); Janicka and Peters (1982); Peters (1984); Frankel et al. (1992b); 
Gaffney et al. (1992); for recent reviews, see Givi (1989); Priddin (1991)). Due to their 
special mathematical properties, the Beta and/or the Dirichlet frequencies yield relatively 
simple analytical solutions for the higher order statistics of the reacting fields. From 
this point of view, the use of the PF is more practical than the AMC since the solution 
procedure does not require the numerical treatment of the PDF transport equation. This 
point has been discussed in detail by Girimaji (1991b). However, as indicated above, the 
use of the Dirichlet frequency cannot be justified for modeling of unpremixed reacting 
flow in a general sense. Also, there is no way of implementing this density directly for 
modeling of non-equilibrium flames, involving strong correlation of the temperature and 
the species mass fractions. Even with the assumption of statistical independence of the 
reacting species and the temperature, the question of the local scalar range imposes a 
severe restriction on the validity of this approximation. For example, it is demonstrated 
by Gaffney et al. (1992) that in the modeling of a reacting turbulent shear flow, depending 
on the a priori choice of the magnitudes of the local scalar bounds the predicted results 
can be altered significantly. Obviously, this problem is not eliminated with the usage of 
JET frequencies in either a univariate or multivariate sense. 

9 CONCLUDING REMARKS 

It is shown that the family of frequencies generated by the Johnson-Edgeworth Translation 
(JET) provides a reasonable means for statistical modeling of binary symmetric scalar 
mixing in homogeneous turbulence. It is also shown that the results predicted by one 
of the members of this family is identical to the solution generated by the Amplitude 
Mapping Closure (AMC) of Kraichnan. This similarity is useful in two regards: (1) 
establishing a mathematical reasoning for the similarity of the probability frequency of 
the Pearson Family (PF) and that of the AMC for the description of the problem, and (2) 
suggesting the possible use of other members of the JET frequencies in approaches in 
which the Probability Density Function (PDF) is assumed a priori The PDF’s generated 
by all these models are shown to compare well with each other and also with the results 
obtained by Direct Numerical Simulations (DNS). However, none of the models are 
capable of accurately predicting the conditional expected dissipation and the conditional 
expected diffusion of the scalar field. This problem is associated with the incapability of 
the models to account for the migration of the scalar bounds as mixing proceeds. 
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